AD-Alll  872 

UNCLASSIFIED 


OHIO  STATE  UNI V  .COLUMBUS  DEPT  OF  GEODETIC  SCIENCE  A— ETC  F/G  8/5 
GLOBAL  GEOPOTENTIAL  MODELLING  FROM  SATELLITE-TO-SATELLITE  TRACK— ETC  (U> 
OCT  81  0  L  COLOMBO  P19628-79-C-0027 

317  AFGL-TR-81-0319  NL 


DIK  FILE.  COEy  AD  A  1 1 1 8  7  2 


AFGL-TR-81-0319 


GLOBAL  GEOPOTENTIAL  MODELLING  FROM  SATELLITE -TO-SATELLITE  TRACKING 


Oscar  L.  Colombo 


Department  of  Geodetic  Science  & 
Surveying  ^ 

The  Ohio  State  University 
1958  Neil  Avenue 
Columbus,  Ohio  43210 


October,  1981 
Scientific  Report  No.  10 

Approved  for  public  release;  distribution  unlimited 


AIR  FORCE  GEOPHYSICS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
UNITED  STATES  AIR  FORCE 
HANSCOM  AFB,  MASSACHUSETTS  01731 


0O  03  09  021 


UNCLASSIFIED  _ 


SECURITY  CLASSIFICATION  OF  this  rage  flWi#/i  0«lEtil>n« 


REPORT  DOCUMENTATION  PAGE 


TO 


i.  report  number 

AFGL-TR-81-0319 


4.  TITLE  (Mtd  Submit) 

GLOBAL  GEOPOTENTIAL  MODELLING  FROM  SATELLITE-TO- 
SATELLITE  TRACKING 


7.  AU  THORf t) 


Oscar  L.  Colombo 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


3-  RECIPIENT'S  CATALOG  NUMBER 


5.  TYPE  OP  REPORT  &  PERIOD  COVERED 

Scientific  Report  No. 10 

6~  PERFORMING^)  1C.  REPORT  NUMBER 

Dept.ofGeod.  Sci.No.317 


6.  CONTRACT  OR  grant  N U M 8 E Rf 47 


F1962G- 79-C-0027 


9.  PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

Department  of  Geodetic  Science  &  Surveying 
The  Ohio  State  University  -  1958  Neil  Avenue 
Columbus,  Ohio  43210 


II.  CONTROLLING  OFFICE  NAME  ANO  AOORESS 

Air  Force  Geophysics  Laboratory 

Hanscom  AFB,  Massachusetts  01730 

Contract  Monitor :  George  Hadgigeorge  /LWG 


4.  MONITORING  AGENCY  NAME  A  ADDRESS/1//  di/larmnt  from  Controlling  Ollica)  15.  SECURITY  CLASS,  (of  thla  raport) 

Unclassified 


10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  4  WORK  UNIT  NUMBERS 

611J2F 

2309G1AW 


12.  REPORT  DATE 

October  1981 


13.  NUMBER  OP  PAGES 

137 


1«.  DISTRIBUTION  STATEMENT  (of  thla  Raport) 


15#.  DECLASSIFICATION/ DOWNGRADING 
SCHEDULE 


A  -  Approved  for  public  release;  distribution  unlimited 


17.  DISTRIBUTION  STATEMENT  (of  tho  abstract  an  farad  In  Block  2D,  II  dillorant  Irons  Raport) 


IS.  KEY  WORDS  ( Contlnua  on  raver  me  aido  If  nacaaaary  and  Idanllly  by  black  number) 


geodesy,  gravity,  satellite-to-satellite  tracking 


20.  ABSTRACT  (Continue  on  raver aa  aide  U  nacaaaary  and  Idanllly  by  black  number) 

'l( The  error  analysis  of  the  global  modelling  of  the  geopotential  has  been 
carried  out  up  to  degree  and  order  331  of  the  spherical  harmonic  expansion, 
for  data  from  a  low-low  satellite-to-satellite  tracking  (SST)  mission.  The 
sphericity  and  the  rotation  of  the  Earth  have  been  considered,  as  well  as  the 
discrete  nature  of  the  data,  assumed  to  consist  of  time  averages  of  the  measured 
range-rate  sampled  at  regular  intervals.  The  expansion  of  the  potential  has  been 
truncated  at  degree  n  *  331,  because  little  information  on  higher  degrees  is 
likely  to  be  present  in  the  data.  Two  theories  have  been  used:  that  of  least  ^ 


DD  i  ISU’L  1473  EDITION  OF  1  NOV  «S  IS  OBSOLETE 
*'  1  JAN  71  nvrt  AOOTDTPn 


EDITION  OF  I  NOV  «S  IS  OBSOLETE 


SECUNITY  CLASSIFICATION  OF  THIS  PAOS  (***>  Dmlt  BnltrtB) 


_ UNCLASSIFIED _ 

stcuwtTv  classification  of  this  PAoen m«n  o«<«  Cntmt»a) 

*  Squares  adjustment,  and  that  of  least  squares  collocation;  above  degree  n  =  zuu 
the  accuracies  predicted  according  to  collocation  are  significantly  better  than 
those  according  to  least  squares.  In  this  report  there  is  also  a  discussion 
on  how  to  process  SST  data  to  obtain  very  high  resolution  models  of  the  grav¬ 
itational  field.  Descriptions  and  listings  of  computer  programs  are  included. 

To  reduce  the  computer  time  and  storage  needed  to  set  up  and  to  invert  the 
normal  matrix,  a  somewhat  simplified  orbital  geometry  and  an  approximate  model 
of  the  data  have  been  adopted;  no  orbit  determination  errors  have  been  considere 
Some  arguments  are  given  to  justify  these  shortcuts,  which  may  not  affect  ser¬ 
iously  the  validity  of  the  results.  An  extention  of  the  theory  to  non-polar 
orbits  is  given. 

The  main  results  according  to  collocation  can  be  summed  up  as  follows:  if 
the  two  satellites  move  in  much  the  same  polar,  circular  orbit  at  a  height  of 
160  km  and  at  a  distance  of  300  km  from  each  other;  if  the  accuracy  of  the  aver¬ 
aged  range  rate  is  /?  x  10”6m  s-1  ,  the  averaging  interval  is  4  s  ,  and  sampling 
takes  place  every  4s;  if  residual  data  are  used,  with  respect  to  a  reference 
model  of  specified  accuracy,  complete  to  degree  and  order  20,  then: 

(1)  the  relative  error  in  the  estimated  potential  coefficients  could  be 
better  than  1%  up  to  degree  n  *  130  ,  than  10%  up  to  n  -  210  ,  and  than  50% 
up  to  n  *  270  ; 

(2)  the  accuracy  of  point  geoid  undulation  implied  by  the  coefficients 
could  be  better  than  0.05  mm  rms  in  the  band  from  3000  km  to  40030  km  (total 
error  in  this  band),  and  better  than  10  cm  rms  in  the  band  from  140  km  to  3000  km 
(also  total  error). 


JCCUWTY  CLASSIFICATION  OF  »</'■  AOCPwim  0.(.  Er 


Foreword 


This  report  was  prepared  by  Or.  Oscar  L.  Colombo.  Much  of  the  research 
described  in  this  report  was  carried  out  while  Dr.  Colombo  was  a  Post 
Doctoral  Research  at  The  Ohio  State  University  where  the  studies  were 
supported  under  Air  Force  Contract  No.  F19628-79-C-0027,  The  Ohio  State 
University  Research  Foundation  Project  711664.  The  contract  covering 
this  research  is  administered  by  the  Air  Force  Geophysics  Laboratory, 

Hanscom  Air  Force  Base,  Massachusetts,  with  Mr.  George  Hadgigeorge,  Contract 
Monitor.  The  actual  writing  of  this  report  was  carried  out  by  Dr.  Colombo  . 
at  the  Geodetic  Institute  of  the  University  of  Stuttgart  in  the  Federal 
Republic  of  Germany. 


iii 


Accession  For 

"nTIS  G^AJcI 
DTIC  Tf.B 

Unannounced 

Justification 


By - 

Distributer./ 
AvallaV.i ' '  v 

'  a.-  . '"'/or 

!  D  \ t  .' '.’-■•rt’i 


no 


Acknowledgements 


I  wish  to  thank  Professor  Rapp,  Project  Supervisor,  for  providing 
the  encouragement  and  the  support  that  allowed  me  to  begin  and  to  carry 
out  this  research.  Thanks  also  to  Reiner  Rummel ,  who  is  partly  respon¬ 
sible  for  my  interest  in  this  problem,  and  who  provided  useful  comments 
on  parts  of  the  manuscript;  to  Carl  Wagner,  for  some  interesting  corres¬ 
pondence  on  SST;  and  to  Susan  Carroll,  who  did  the  typing. 

This  report  was  written  during  the  first  months  of  a  visit  to  the 
Geodetic  Institute  of  the  University  of  Stuttgart,  in  the  Federal  Republic 
of  Germany.  This  visit  was  sponsored  by  the  Alexander  von  Humboldt  Founda¬ 
tion.  To  the  staff  of  this  institution,  and  to  Professor  Erik  Grafarend, 
whose  liberality  allowed  me  to  finish  this  work,  my  most  sincere  appreciation. 


v 


Table  of  Contents 


Foreword . iii 

Acknowledgements  . 

1.  Introduction  .  1 

1.1  The  Low-Low  Configuration  .  5 

1.2  The  Band  Limited  Assumption . 6 

1.3  An  Approximate  Model  for  the  Line  of  Sight  Velocity  .  12 

2.  The  Mathematical  Model . 17 

2.1  Simpl ifying  Assumptions  .  17 

2.2  The  Extended  Legendre  Function . 18 

2.3  Time  Series  Expression  of  the  Inertial  Line  of  Sight 

Acceleration . 20 

2.4  The  Correction  Term . 22 

2.5  The  Observation  Equation . 22 

2.6  The  Condition  that  Nr  and  Nq  Must  Be  Relative  Primes.  .  .  25 

(a)  Case  where  m  f  q . 25 

(b)  Case  where  m  =  q . 25 

2.7  The  Scalar  Product  of  Two  Columns  of  the  A  Matrix . 26 

2.8  Least  Squares  Adjustment . 28 

2.9  The  Structure  of  the  Normal  Matrix . 31 

2.10  The  Existence  of  G"1 . 33 

2.11  Least  Squares  Collocation . 35 

2.12  Accuracy  of  the  Computed  Geoidal  Heights  .  39 

2.13  The  Effect  of  Some  Mission  Parameters  on  Coefficient 

Accuracy . 42 

2.14  The  Right  Hand  Sides  of  the  Normals . 44 

2.15  Oblique  Orbits . 46 

3.  Numerical  Results . 50 

3.1  Spectral  Model  and  Error  Formulas  .  50 

3.2  Results  According  to  Least  Squares  Adjustment  .  52 

3.3  Results  According  to  Least  Squares  Collocation . 54 

3.4  Accuracies  of  Different  Harmonics  of  the  Same  Degree . 63 

4.  Validity  of  the  Results . 64 

4.1  The  Geometry  of  the  Real  Orbit . 64 

4.2  Vertical  Reduction  to  the  Mean  Sphere  .  66 

4.3  The  Effect  of  Errors  in  the  Calculated  Orbits . 73 

5.  Data  Processing . 74 

5.1  An  Iterative  Approach  .  74 

5.2  Other  Methods  .  77 

5.3  The  Use  of  Local  Solutions . 78 

6.  Conclusions . 78 

vii 

raicsfiun  mx  UAM-mr  nitm 


References  .  81 

Appendix  A  :  Orbital  Perturbations  .  34 

Appendix  B:  Computer  Programs  .  88 

B.l  Main  Program . 88 

(a)  Full  Version . 88 

(b)  Reduced  Version . 92 

8.2  Subroutine  ONEREV  .  92 

B.3  Subroutines  LEGFDN ,  MODEL,  and  NVAR . 93 

B.4  Sample  Output  .  94 

Appendix  C:  Detailed  Listings  Degree  by  Degree . 123 


vlii 


1.  Introduction 


More  than  a  decade  ago,  in  1968,  Muller  and  Sjogren  published  their 
paper  on  lunar  mass  concentrations,  or  "mascons".  In  it  they  announced 
one  of  the  most  important  conclusions  about  the  internal  structure  of  a 
member  of  the  Solar  System  ever  reached  from  the  analysis  of  gravitation 
alone.  With  remarkable  insight,  they  used  the  numerically  differentiated 
range-rate  tracking  data  of  the  lunar  orbiters,  deployed  as  part  of  the 
Apollo  program,  as  "direct  observations"  of  the  acceleration  of  gravity 
along  the  line  of  sight  between  the  tracking  station  on  Earth  and  the 
spacecrafts  circling  the  Moon.  A  simple  plot  of  this  information  showed 
strong  and  roughly  circular  anomalies  over  the  flatlands,  or  maria,  sug¬ 
gesting  the  presence  of  large  bodies  of  abnormal  density  buried  in  the 
lunar  crust. 

The  idea  that  high  resolution  information  on  a  gravitational  field 
could  be  extracted  by  a  simple  analysis  from  differentiated  range-rate 
data  was  seen,  naturally  enough,  as  having  potential  value  for  the  study 
of  our  own  planet.  One  more  or  less  obvious  adaptation  of  the  idea  was 
having  two  satellites:  one  in  a  very  high  orbit  tracking  the  second  one 
on  a  low  orbit,  much  as  the  "high"  station  on  Earth  had  tracked  the  "low" 
orbiters  near  the  Moon.  A  different  satellite-satellite  tracking  (SST) 
configuration,  later  to  be  known  as  the  "low-low"  approach  (the  first  being 
the  "high-low",  of  course),  was  proposed  by  Wolff  in  1969:  he  thought 
that  the  relative  velocity  along  the  line  of  sight  between  two  bodies  fol¬ 
lowing  as  close  as  possible  the  same  near  circular,  polar  orbit  was, as 
a  first  approximation,  the  difference  in  anomalous  potential  between  the 
two.  As  the  Earth  rotates,  this  pair  would  have  eventually  covered  the 
whole  planet  with  observations  of  gravity,  all  taken  with  the  same 
"instrument",  creating  a  global  data  set  of  extraordinary  density  and  of 
uniform  quality.  From  such  data  one  should  have  been  able  to  recover  an 
almost  as  detailed,  and  far  more  complete,  picture  of  the  field  than  from 
terrestrial  data  alone.  Though  this  idea  suffered  from  some  theoretical 
and  practical  problems,  its  appeal  to  the  imagination  was  such  that  it 
inspired  much  research  in  spite  of  its  shortcomings. 

Among  the  first  to  help  clarify  the  theory  behind  SST  was  Schwarz 
(1970),  who  proposed  a  rigorous  mathematical  formulation.  Many  studies 
followed  to  find  out  the  best  satellite  arrangement,  how  the  data  could 
be  processed,  and  what  accuracy  could  be  expected  from  the  results.  To 
all  this  one  must  add  the  work  done  on  the  development  of  suitable  hard¬ 
ware,  also  along  several  lines. 

In  the  mi d- 1970' s,  during  the  Apollo-Soyusz  and  the  Geos-3  missions, 
"high-lov."  range-rate  tracking  data  were  gathered  using  the  ATS-6  satellite, 
in  geostationary  orbit,  as  the  high  spacecraft,  because  of  its  steerable 
radar  antenna.  During  the  same  period,  the  idea  of  a  dedicated  gravitational 
satellite  mission  begun  to  develop.  Both  NASA  and  its  European  counterpart 
ESA  drew  up  preliminary  plans  for  a  Gravity  Satellite  (GRAVSAT) ,  and  for 
a  Space  Laser  Low  Orbit  Mission  (SLALOM),  respectively.  Of  the  two,  SLALOM 
is  the  most  concretely  defined  at  present  (CSTG  Bulletin  No. 2,  1980).  It 
involves  the  simultaneous  tracking  by  laser  interferometry  from  the  Space 
Laboratory, carried  on  the  Space  Shuttle,  of  two  reflecting  spheres  to 
determine  their  relative  velocities  with  respect  to  the  shuttle  and  to 
each  other.  Recovery  of  gravity  information  is  to  be  limited  to  a  region 
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over  the  Eastern  Mediterranean,  and  to  a  period  of  some  seven  days.  This 
experiment  is  expected  to  be  carried  out  during  this  decade.  GRAVSAT, 
on  the  other  hand,  is  a  global  concept  that  may  take  one  of  several  possible 
configurations  based  on  the  SST  principle  or,  alternatively,  resort  to 
an  orbiting  gradiometer.  At  present  it  appears  likely  that  the  SST  idea 
will  be  tried  first,  in  a  GRAVSAT-A  mission  in  the  late  1980's,  while  the 
gradiometer,  whose  development  into  a  practical  instrument  for  this  purpose 
is  still  at  the  "breadboard"  stage,  may  be  used  in  a  GRAVSAT-B  mission 
sometime  in  the  1990's.  Several  types  of  gradiometers  have  been  considered. 

A  promising  design  seems  to  be  a  supercooled  instrument  now  under  development 
at  the  University  of  Maryland,  which  may  achieve  a  sensitivity  of  better 
than  0.001  E  (Paik,  1981). 

Of  the  various  implementations  of  the  "low-low"  principle,  the  most 
likely  to  be  adopted  appears  to  be  the  DISCOS  system  of  two  drag-compensated 
satellites,  capable  of  remaining  in  orbit  for  up  to  six  months  so  close 
to  the  Earth  that,  without  the  periodical  use  of  small  rocket  engines, 
their  orbits  would  decay  very  quickly  due  to  the  atmospheric  resistance. 

Two  proof-masses,  one  inside  each  craft,  will  be  kept  in  permanent  free 
fall  by  the  compensating  mechanism,  so  only  gravitational  forces  act  upon 
them.  Their  relative  line  of  sight  velocity  will  be  measured  by  an  extremely 
accurate  radar  interferometric  technique  developed  at  the  Applied  Physics 
Laboratory  of  Johns  Hopkins  University  (Pisacane  and  Yionoulis,  1980). 
Accuracies  of  better  than  10“6ms_1  are  expected,  provided  that  no  serious 
problems  due  to  ionospheric  propagation  are  encountered. 

Apart  from  the  "high-low"  and  the  "low-low"  configurations,  an  intermed¬ 
iate  “butterfly"  arrangement,  where  two  satellites  follow  different  ellip¬ 
tical  orbits,  has  been  considered  as  well. 

The  SST  data  collected  during  the  middle  of  the  last  decade  were  analysed 

in  different  ways.  Among  others,  Kahn  et  al.  (1978)  tried  to  recover  gravity 

anomalies  using  the  "numerical"  method  of  satellite  geodesy,  with  the  aid 
of  the  computer  programme  GEOOYN;  Hajela  (1978)  followed,  instead,  the 
original  idea  of  treating  the  differentiated  range-rates  as  gravity  obser¬ 
vations,  and  used  least  squares  collocation  as  the  processing  technique. 

Marsh  and  Marsh  (1978)  attempted  what  was,  in  essence,  an  experiment  like 
that  of  Muller  and  Sjogren  with  Earth  data,  to  detect  crustal  and  upper 
mantle  structures.  All  these  studies  have  been  restricted  to  local  areas, 
as  no  complete  global  set  of  SST  has  been  obtained  yet. 

The  work  that  has  been  done  on  the  error  analysis  of  the  recovery 
of  values  of  mean  gravity  anomalies,  mean  and  point  geoidal  undulations, 

etc.,  can  be  divided,  broadly,  along  two  main  lines:  the  "numerical"  method 

of  satellite  geodesy,  and  the  "direct"  approach  that  goes  back  to  Muller 
and  Sjogren:  using  the  differentiated  range-rate  values  as  observations. 

As  an  example  of  the  first  line,  one  can  mention  the  work  of  Douglas  et 
al.  (1980),  as  one  of  the  most  recent.  Of  the  second,  the  author  is  best 
familiar  with  the  work  of  Hajela  (1974),  Rummel  et  al.  (1976),  Krynski 
(1978),  Rapp  and  Hajela  (1979),  and  Rummel  (1980).  All  of  these  have  consid¬ 
ered  local  gravity  field  model  improvements  relying  on  the  theory  of  least 
squares  collocation.  The  ultimate  accuracy  of  global  recovery  given  a  certain 
quality  of  data  has  been  investigated,  among  others,  by  Breakwell  (1979), 
who  used  a  flat-earth  approximation,  while  Jekeli  and  Rapp  (1980)  have 
employed  a  spherical,  non-rotating  Earth  approximation. 
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While  all  the  work  mentioned  above  has  been  in  progress,  numerous 
meetings,  involving  government  agencies  and  members  of  the  scientific  com¬ 
munity  that  are  potential  users  of  GRAVSAT  data,  have  taken  place.  Of 
several  reports  on  these  activities,  there  are  two  of  particular  importance 
that  describe  the  objectives  of  a  GRAVSAT  mission  and  also  define  the 
required  accuracies  and  other  specifications  for  the  results  so  they 
can  be  used  meaningfully  for  geodetic  and  geophysical  purposes.  One  is 
the  report  of  a  special  panel  of  the  Committee  on  Geodesy  of  the  US  National 
Academy  of  Sciences  (1979);  the  other  is  by  the  GRAVSAT  Users  Working 
Group  (1980).  The  main  requirements  defined  so  far  are: 

a)  For  geological  and  geophysical  applications,  qrav’ cy  information  should 
be  resolved  at  the  Earth's  surface  through  wavelengths  of  100  km  and  to 

an  accuracy  of  2.5  to  10  mgals; 

b)  Oceanographers  need  geoid  heights  accurate  to  some  10  cm  in  the  band 
from  100  km  to  3000  km. 

The  first  specification  is  relevant  to  the  study  of  the  structure 
of  the  crust  and  the  upper  mantle.  The  second  one  is  associated  with  the 
analysis  of  the  instantaneous  and  average  shape  of  the  sea  surface  in  studies 
that  use  such  data  as  satellite  altimetry.  The  difference  between  the 
sea  surface  and  a  horizontal  surface  can  reveal  many  as  yet  unknown  aspects 
of  surface  currents,  tides,  and  transients  of  various  kinds,  particularly 
in  regions  of  the  oceans  that  are  not  sufficiently  accessible  by  other 
means.  Also  of  consequence  to  oceanography,  as  well  as  to  many  practical 
aspects  of  geodesy,  such  as  satellite  positioning  techniques,  is  the  obtention 
of  a  model  for  the  gravity  field  that  permits  the  calculation  of  very  precise 
spacecraft  orbits,  so  this  is  another  important  objective  of  GRAVSAT.  Among 
the  latest  global  studies  those  by  Breakwell ,  and  by  Jekeli  and  Rapp, 
already  mentioned,  suggest  that  the  DISCOS  system  could  achieve  both  the 
quality  and  quantity  of  data  needed  to  meet  goals  (a)  and  (b). 

The  present  report  describes  the  theory  behind  a  global  error  analysis 
of  a  low-low  mission  of  the  DISCOS  type,  and  gives  the  corresponding  results. 
The  latter  are  in  broad  agreement  with  some  of  the  previous  studies,  notably 
Breakwell's  and  Jekeli  and  Rapp's.  The  approach  taken  has  been  the  "direct" 
one  of  considering  the  differentiated  range-rate  signal  as  equal  to  the 
component  of  the  gravitational  line  of  sight  acceleration,  in  an  inertial 
frame,  along  the  line  of  sight  direction.  This  is  only  an  approximation, 
but  it  simplifies  greatly  the  mathematical  treatment,  and  previous  reports 
by  Hajela  (1978)  and  by  Rummel  (1980), already  mentioned,  indicate  that 
it  may  be  in  good  agreement  with  reality  at  short  spatial  wavelengths, 
both  for  the  "high-low"  and  the  "low-low"  configurations,  respectively. 

This  model  of  the  line  of  sight  signal  has  been  modified  somewhat  here 
by  substracting  a  term  that  depends  only  on  radial  distance  to  the  geocenter 
to  eliminate  some  unrealistic  long-wave  phenomena  related  to  the  large 
zero  harmonic  and  to  the  other  even  zonal s  of  the  geopotential. 

This  work  is  concerned  primarily  with  the  optimal  recovery  of  a  geo¬ 
potential  model  in  the  form  of  a  spherical  harmonic  expansion  of  the  potential 
of  such  a  high  degree  and  order  (331)  that  its  truncation  error 
at  satellite  altitude  (160  km)  is  nearly  negligible,  particularly  in  the 
presence  of  noise. x1)  It  differs  from  previous  studies  of  this  sort  in 
that  it  considers  a  spherical,  rotating  Earth,  and  discrete  data  consisting 
of  time-averages  of  the  instaneous  range  rate.  Observation  equations  are 
obtained  under  the  simplifying  assumptions  that  the  orbits  are  perfectly 
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circular,  the  satellites  separation  is  constant,  and  that  the  orbits  repeat 
themselves  exactly  every  179  days,  i.e.,  the  length  of  the  whole  mission. 

From  these  equations  a  normal  matrix  is  formed  that  is  then  inverted 
taking  advantage  of  its  block-diagonal  structure.  Two  adjustment  techniques 
are  considered  and  implemented:  least  squares  adjustment,  and  least  squares 
collocation.  An  argument  based  on  an  infinite  series  of  successive  approxi¬ 
mations  is  used  to  show  how  the  results  obtained  under  such  simplifying 
assumptions  can  be,  in  fact,  close  indicators  of  the  maximum  amount  of 
information  than  can  be  extracted  from  a  real,  three-dimensional  distribution 
of  SST  data.  While  orbital  errors  are  not  included  in  the  analysis,  it 
turns  out  that  the  effect  of  such  errors  on  the  estimated  parameters  is 
very  small,  although  this  particular  result  depends  strongly  on  the  model 
adopted  for  the  line  of  sight  signal. 

The  last  section  considers  how  the  principles  used  in  this  study  can 
be  applied  to  the  actual  processing  of  SST  data  in  order  to  obtain  a  very 
high  resolution  spherical  harmonic  model  of  the  gravity  field.  Such  models 
have  important  mathematical  advantages,  and  there  is  no  great  problem  in 
using  them  if  this  is  done  with  adequate  techniques  (see  Colombo 
(1981),  for  instance).  Global  data  reduction  is  a  very  important  problem, 
to  which  not  enough  attention  has  been  paid  so  far. 

The  potential  impact  of  SST  and  of  satellite  gradiometry  on  the  future 
of  geodesy  and  of  geophysics  appears  great.  It  is  a  difficult 
and  beautiful  task,  the  work  of  many,  through  many  years,  to  develop 
what  begun  as  a  simple  and  bright  idea  for  studying  the  Moon  into  a  tool 
for  increasing  our  understanding  of  our  own  planet  and,  eventually, 
of  the  rest  of  the  Solar  System.  If  we  are  successful,  this  task  could 
have  a  deeper  and  more  lasting  effect  on  pure  and  applied  Earth  sciences 
than  any  other  ever  attempted  in  geodesy  before. 


In  this  "band-limited"  situation,  the  optimal  estimates  of  any  other 
field-related  quantities  (undulations,  gravity  anomalies,  etc.)  can  be 
obtained  directly  from  the  optimal  estimates  of  the  potential  coefficients 
(see,  for  instance,  Colombo  (1981),  paragraph  2.18).  The  same  is  true  of 
the  accuracies  of  those  quantities. 
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1 . 1  The  Low-low  Configuration 

Figure  1.1  shows  two  drag-compensated  satellites  S:  and  S2  that  have 
been  placed  almost  in  the  same  near-circular  polar  orbit  so  that,  in  inertial 
space,  the  plane  of  the  common  orbit  contains  the  Earth's  mean  axis  of 
rotation,  which  is  aligned  in  the  picture  with  the  xa  axis.  The  geocentric 
angle  separating  both  satellites  is  &  =  2  sin"^-^-)  radians,  where  0 
is  the  length  of  the  segment  of  1  ine  of  sight  between  b :  and  S?  .  The  positions 
of  the  spacecrafts  in  the  system  of  geocentric  inertial  crordinates  xi, 
x.,  Xj  are  represented  by  the  vectors  x:  and  x,  ,  respectively.  The  line 
of  sight  vector  Xi;  =  xi-xi  is  oriented,  according  tc  the  picture,  from 
South  to  North  in  the  ascending  passes.  The  sense  of  the  orbit  (prograde 
or  retrograde)  is  not  important.  The  mean  orbital  radius  is  R  . 

Both  satellites  turn  round  the  Earth  with  approximately  the  same  angular 
velocity  w  =  -/GMR' ?  while  the  Earth  itself  turns  on  its  axis  with  mean 
angular  velocity  Q  .  The  points  directly  below  each  satellite  describe 
groundtracks  that  envelope  the  whole  surface  of  the  planet  as  the  mission 
progresses.  Both  satellites  have  the  same  instantaneous  longitude,  but 
their  groundtracks  are  not  identical:  they  are  always  '■P rad  apart  in  the 
S-N  direction  and  tprad  in  the  E-W  direction.  In  what  follows,  the 
word  "groundtrack" ,  unless  further  explanation  is  given,  should  refer 
to  that  of  the  point  midways  between  the  satellites,  along  the  line  of  sight, 

P  “  i  ill 2  • 

The  measuring  device  detects  the  relative  velocity  along  the  line 
of  sight  between  two  proof -masses ,  one  inside  each  satellite,  which  are 
kept  in  permanent  free-fall  round  the  Earth  by  the  drag-compensating  mechanism. 
This  mechanism  compensates  not  only  for  drag,  but  for  all  other  non-gravi- 
tational  influences  of  consequence  as  well.  The  relative  line  of  sight 
velocity  is 

V12  =  Al2  !l2  (I-l) 

where 

A12  =  Ai  "  *2 

is  the  relative  inertial  velocity,  x,  and  x0  beinq  the  absolute  ones, 
while  1 

6.1 2  ~~P  1  ill 2  (1-2) 

is  the  unit  vector  along  the  line  between  and  Sj  ,  and 

P=  lliisll  =  /(XU-XU)'J  +  (X21-x22)2  +  (X31-X32/'2_  d-3) 

is  the  euclidean  norm  of  Xia  >  i.e.,  the  distance  between  the  two  spacecrafts. 
The  observed  values  of  v12  are  averaged  over  Aa  seconds  and  then  transmitted 
to  earth-side  stations  every  At  seconds,  where  Aa  <  At  .  These  averages 
constitute  the  SST  data  to  be  studied  hr -e,  and  the  expression  for  any 
one  of  them  is 

v12(t)  =  Aa  Vi2(t)  dx  (1.4) 
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The  line  of  sight  velocity  v12  reflects  small  changes  in  the  velocities 
of  both  satellites  about  their  common  average  u  =  Rw  .  These  changes 
are  brought  about  by  the  anomalous  gravitational  field,  which  is  the  dif¬ 
ference  between  the  true  field  and  that  of  a  homogeneous  geocentric  sphere 
of  the  same  mass  as  the  Earth  and  radius  smaller  than  R  . 


t 

L 


The  actual  gravitational  potential  V  can  be  represented  in  geocentric 
spherical  coordinates  (radial  distance  r  ,  latitude  $  ,  longitude  X  ) 
by  a  spherical  harmonic  expansion 


V(r,  p,  X)  =  Jq  Pnm(sin*)  cosmX+  5nm  sinmX]  (1.5. a) 

where  P  (sirvj))  is  the  fully  normalized  associated  Legendre  function  of 
the  first  kind,  and  Cp,m  5nm  are  fully  normalized  spherical  harmonic 
coefficients.  The  following  alternative  notation  will  be  used,  wherever 
possible,  in  this  work: 

«'•♦•»>■  *  1 1 <>•»> 
n*u  m*u  a=0 

Wlth  Ca  7°  (<t>,X)  -  P  (sino)  |  Snm  c(?sm*  a=? 
nm  nm  '  ’  '  nm  '  \  5_m  sinmX  a  =  1 

nm 


M  is  the  mass  of  the  Earth,  6  the  universal  constant  of  gravitation, 
"a"  is  the  mean  Earth  radius,  "n"  indicates  the  degree,  and  "m"  the 
order  of  each  term  in  the  expansion. 
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The  three  terms  with  n=l  are  zero,  because  here  the  geocenter  coincides 
with  the  origin  on  coordinates;  the  zero  harmonic 
equals  ,  so  Coo=  1  .  The  anomalous  potential  is,  therefore, 

R(r,*,A)  =  V(n»,A)  -  (1.6) 


The  disturbing  potential  T  ,  the  modelling  of  which  is  the  concern  of 
this  report,  is  the  difference  between  the  true  potential  V  and  some 
reference  model  potential  U  of  the  form 

NM  n  1 


i.  _  GM 

U  (r,i $»a)  -  — 


l  l  l 

n=0  m=0  a=0 


■  a  (M) 


nm 


(1.7) 


where  NM  is  a  relatively  small  integer  (20  or  30).  The  objective  of 
this  study  is  finding  out  the  accuracy  with  which  a  model  of  T  of  the 
same  form  as  U  ,  but  truncated  at  a  much  higher  order  N  ,  could  be  recovered 
from  low- low  SST  data. 


As  explained  later  in  this  section,  the  time  derivative  of  the  line 
of  sight  velocity  can  be  approximated  by  the  component  of  the  inertial 
acceleration  aimed  along  the  line  of  sight,  minus  a  term  independent  from 
0  and  \  .  This  acceleration  is  a  linear  combination  of  the  three  accel¬ 
erations 


ar  (r,0 ,x) 
%  ( r,<M ) 
a  x  (r,(j> ,x ) 


=  11  =  y  y  y  f —f  ca  9°  (p.x) 

3r  nW0a=0  “  nm  nm 

_  1  3V  .  GM  r  r  p  (jLf1  ha  JL-  7a 

r  3«  r‘  n=0  m=0  a=0  ^  nm  3*  nm 

1  ll  =  -  r  r  r  fJLf  Ca  ya  (1  8,c) 

3  A  r*c5 s<*  nL  mL  L  nm  9  a  nm 


(1.8,a) 
( 1 .8 ,b) 


r  cos  <p 


where 


J_  7a  =  J-  P 
3d  nm  a  a  nm 


mx 


and 


i-  r  =  m  P 

3A  nm  nm 


m\ 


3ecause  all  these  expansions  converge  outside  the  Earth’s  bounding  sphere, 
the  general  terms  of  all  of  them  should  tend  to  zero  with  n  tending  to 
infinity,  so  the  size  of  the  harmonics  should,  in  general,  decrease  as 
n  increases.  This  decay  must  be  accentuated  by  the  factors  (f-)n  ,  where¬ 
fore  the  field  should  become  smoother  with  altitude,  as  the  higher  frequency 
terms  vanish  faster  than  the  rest  with  increasing  r  .  At  satellite  height 
h=R-a  this  smoothing  should  mean  that,  above  a  certain  degree  N  ,  all 
terms  in  the  expansion  can  be  neglected.  Consequently,  the  field  can  be 
regarded  as  band- 1 imi ted,  with  terms  restricted  to  degrees  in  the  band 
0  *  n  c  N  .  As  the  error  analysis  method  presented  in  section  2  requires 
computer  time  in  proportion  to  N1'  ,  it  is  important  to  find  a  realistic 
value  for  N  that  is  also  as  low  as  possible.  The  reasoning  that  follows 
attempts  to  provide  a  guide  for  such  a  choice  using  the  decay  in  the 
spectrum  of  the  line  of  sight  inertial  acceleration  as  a  criterion. 


To  simplify  matters.  Earth  rotation  can  be  ignored,  the  orbit  can  be  assumed 
to  be  perfectly  circular,  and  the  geocentric  separation  i>  to  be  constant. 
Under  these  conditions,  all  field  dependent  functions  are  periodical,  the 
line  of  sight  acceleration  among  them,  and  can  be  represented  by  Fourier 
series  such  as  , 


a 1 2 (t )  =  l  a.  coskajt  +  bk  sinku>t 
k=0 

As  the  Earth  does  not  rotate  in  this  case,  one  can  choose  an  arbitrary 
system  of  geocentric  coordinates  r'=r  ,  ■$>',  X'  where  the  "equator"  coin¬ 
cides  with  the  plane  of  the  orbit,  and  then  make  the  substitution 

X  =  MODULE  2tt 


so  that  the  Fourier  series  becomes 


a i ,{ A)  =  l  at  cos kX  +  bk  sinkX 
k=0 

To  find  out  the  coefficients  of  this  series,  consider  first  the  inertial 
line  of  sight  acceleration  in  this  system  of  coordinates.  For  the  first 
satellite,  let 


/ ar  (R,  <y  =  0,  X')  \T  T 
a j.  =  {  a 4>  (R,  <y  =  0,  X')  )e.i2  =  £i  £12 
VaX  (R,  <y  =0,  X')  / 

be  the  projection  of  its  inertial  acceleration  along  that  line,  and  let 
a  2  =  5.1  ®  1 2 

be  the  corresponding  projection  for  the  second  satellite.  Then  the  line 
of  sight  acceleration  is 


aia^L(S(Ss)rPO+V  (Sx)^> 
-  et,(ar(S2)  rSSa)  +  a^.  (S2)  $iSz) 


+  ax,  (Si)  x[Sl)) 
+  ax,  (S2)  x[Sz)) 


(1.9) 


where  r^,  and  X^  are  the  unit  vectors  pointing  upwards, "S-N" , 

and  "W-E"  at  the  general  point  P(r,<f>'  ,X' ) .  Because  the  circular  orbit  new 
lies  on  the  equatorial  plane  of  the  rotated  system  of  coordinates,  we  have 


el*'5'’  •  =I2*iS‘>  •  0 


(l.lO.a) 


From  Figure  1.2  it  is  easy  to  see  that 


eLrSS‘>  • 


(l.lO.b) 


and 


T  AS,)  ,  .T.,1  S2) 


6 1  2  ^ 


6.1  2Ao 


(l.lO.c) 


Calling  the  longitudes  of  Si  and  S2  X'  and  X'  +  v  ,  respectively, 
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replacing  (1.8,  a-c)  and  (1.10,  a-c)  in  (1.9),  rearranging  terms  (which 
is  valid,  because  all  the  expansions  converge  uniformly  outside  the  bounding 
sphere),  and  making  use  of  simple  trigonometric  identities, 

ai2(R,*’=0,  V)  =  fy  l  I  5nm(0)(^n  C(n+l)Cnmsin  |-+  m$nmcos  |-]  cosnV 
K  in=0  n=m 

+  [(n+l)5nnisin|  -  mCpmcos|- ]  sinmX1  -  [- (n+1  )Cnm s i n  \  +  mSnm  cos  ’^3 

.  cos  m  ( \ '  +\p ) 

-  [-(n+1)  5nmsin^-  -  m^nmc0S7'^  sinm(A’+^)  (1.11) 

Therefore,  the  coefficients  of  the  Fourier  series  for  a 12  are 

1  n=m 

+  m  cos  4j-  sin rmp]  +  5nm[m  cos  (1-cos mip)  +  (n+1)  sin  sinm  iii]}  (1.12, a) 
d  n  ct  co 

bm  =  It  l  Pnrn(0)(7)niCnm[-(n+1)sinf  -  m  cos  |  ( 1-cos  m*)]  + 

K  n=m 

+  5nmC (n+1)  sin  (1+cosm^)  +  m  cos  sin }  ( 1 . 12 ,b) 


GM 

am  "  jr 
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where  "m"  has  replaced  the  original  subscript  "k“  for  obvious  reasons. 

The  time  integral  /  ai2(t)  dt  has  a  term  of  the  form  a0t  so  it  cannot 
be  identical  to  v12  because  the  latter  is  periodical  under  the  current 
assumptions.  This  problem  can  be  solved  by  removing  the  constant  term 
from  the  Fourier  series  of  a22  ,  so  the  time  derivative  of  vi2  would 
be  (approximately) 


ai  2  =  2  ■  a0 

S  Vj  2 


(1.13) 


'12=  1  a12  (R.O.X1)  dV  =  f0 


(ait1)  dt* 


=  y  a  SllLJ^  -  b 
L-  m  m  u) 


cosm  mt 
m  m  w 


+  l  2bl 

m=l 


If  (1.13)  can  be  accepted  as  a  valid  approximation,  then 

4-  a^2  dt  =  i  f  a2  +  b2  (where  T,  =  -4-) 

T  Jo  w  _  i  m  m  5 

s  m=l 

00 

=  y  p 

111 

m=l 


(1.14) 


where  Pm  =  i(a2+b2)  is  the  power  at  the  mth  frequency.  Consider  now 
the  averaging  operator  M  {}  ,  which  will  be  encountered  again  in  Section  2 
in  connection  with  the  use  of  least  squares  collocation  as  an  alternative 
to  the  usual  least  squares  adjustment  for  obtaining  a  field  model.  This 
operator  represents  an  average  over  all  rotations;  here  it  can  be  seen 
as  averaging  over  all  possible  circular  orbits  round  a  non-rotating  Earth. 
As  shown,  for  instance,  in  (Colombo,  1981,  Section  2), 


H  ■  «  <5™>  •  Sn+T 

M  (Cnm  Skq  }  =  0  for  all  integer  n  ,  m  ,  k  ,  and  q 


where 


B‘£n»V  sM<SmnV  =  0 

n 

a2  =  y  C2  + 

un  nm  nm 


for  n  /  k  ,  m  f  q 


(1.15,a) 
(1.15, b) 
(1.15,0 


(1.16) 


is  usually  called  the  "nth  degree  variance"  of  the  potential  coefficients. 
The  average  orbital  power  per  wavelength  is,  then. 


(1.17,1) 


for  a12  ,  and 


s  a  -  ..*1  P 
m  arm7  m 


(1.17,b) 


for  v12  .  Squaring  (1.12,  a-b)  and  modifying  the  resulting  expansion 
according  to  (l.l5,a-c),  (1.17,b)  finally  becomes 
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Sm  =  Prm  U+c°smo)  +  Pfc  (1-cos  m>)  +  2P  f  sinrw 


where 


m  1  rm  Vi  r  tm 

n2  m2  “ 


’  "  p™  ’  tot  I  (l>2" (n+1)!  (1-cos*> p™  <°> 

n=m 

is  the  contribution  from  the  radial  acceleration, 

"t.  ■  TI&&-X  Si  <R->2"  ■*<»««>**>  ^  ») 

the  corresponding  contribution  from  the  along-track  acce. leration ,  and 

Prtm  =  ZR\>*  JL  Zn^T  sin*  Pnm(0) 


(1.18) 
(1.19, a) 

(1.19,b) 

(1.19.C) 


is  the  average  crosspectral  power  of  the  two. 


At  ajieight  R-a  =  160  km,  the  angular  frequency  of  the  satellites 
is  u  =  =  1.196  x  10‘3rad.  s*1.  If  the  separation  between  the  satellites 

is  p  *  300  km,  so  «|>  s  0.04594  rad.  ,  ar.d  if  the  an  are  those  used  in 
the  error  analysis  of  Section  3,  empirically  derived  from  terrestrial 
and  satellite  data,  then  the  Sm  are  as  listed  in  the  table  below: 

Table  1.1 

Spectral  R.M.S.  (S  )*  of  the  Line  of  Sight  Velocity, 


III  . 

and  (f>  )*  (acceleration) 


spatial  frequency  m 
(cycles  per  rev.) 


(cm  s  *) 


(mgal ) 


0 

0. 

— 

1 

.335 

.410 

2 

48.854 

119.813 

3 

.348 

1.279 

4 

.205 

1.005 

5 

.170 

1.044 

10 

.678 

X 

10'1 

.833 

20 

.268 

II 

.656 

30 

.192 

II 

.706 

40 

.130 

II 

.637 

50 

.932 

X 

IO-2 

.572 

100 

.855 

X 

io- 3 

.105 

200 

.394 

X 

10-“ 

.967 

X 

10“ ; 

300 

.'107 

X 

IO- 5 

.394 

X 

io- 3 

400 

.207 

X 

IQ'7 

.101 

X 

10* - 

700 

.817 

X 

io- 11 

.701 

X 

10** 

1000 

.544 

X 

io- lu 

.660 

X 

10*11 

Total  rms  of  signal  above 
m  =  0 


48.86  cm  s"1 
0.06  cm  s_1 
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Hgsl 


To  obtain  the  values  listed  above,  the  series  in  ( 1 . 19 ,a-c )  were  truncated 
at  n=110&  though  much  the  same  results  were  obtained  with  n=500  (up  to m  =300), 
which  suggests  a  strongly  band-limited  nature  for  a12  and  Vi2.The  orepon- 

derance  of  the  term  with  m=2  is  due  to  the  large  second  zonal  is 
related  to  the  Earth  oblatness.  The  signal  in  the  data  consists  of  time 
averages  of  vi2  ,  so  its  spectrum  should  be  smoother  than  the  entries 
in  the  table  suggest.  For  the  purpose  of  this  discussion,  such  a  refine¬ 
ment  is  not  necessary,  and  the  effect  of  time-averaging  will  not 
be  considered  until  Section  2. 

With  At=4  seconds,  some  s  1427  samples  of  vi2  are  taken  during 

each  revolution  of  the  satellites.  Accordingly,  the  "Nyquist  frequency" 
of  the  data  is  Ny=713  .  It  is  clear  from  the  table  that  the  power  above 
this  frequency  is'  negl igible,  so  that  aliasing  problems  related  to  the 
sampling  rate  are  likely  to  be  insignificant.  The  question  of  choosing 
N  ,  the  highest  degree  in  a  band-limited  model  of  the  potential,  so  that 
this  model  can  be  regarded  as  realistic,  is  not  an  easy  one,  except  for 
the  fact  that  N  does  not  have  to  be  larger  than  713.  For  this  study 
the  value  of  N=331  has  been  chosen  purely  on  the  basis  that  this  number 
appeared  to  be,  on  preliminary  estimates,  the  largest  N  compatible  with 
the  computing  resources  available  to  the  author.  The  results  in  Section  3 
show  recovery  errors  of  more  than  80%  of  the  actual  values  for  the  coeffic¬ 
ients  of  harmonic  n=330  and,  for  a  number  of  reasons  discussed  in  Section  4, 
these  estimates  are  rather  optimistic  at  the  upper  end  of  the  spectrum. 

So,  perhaps,  N=331  is  truly  close  to  the  upper  limit  of  resolution  for 
global  estimation  procedures  of  the  kind  considered  here  (least  squares 
and  collocation). 


1.3  An  Approximate  Model  for  the  Line  of  Sight  Velocity 

To  estimate  the  accuracy  with  which  the  spherical  harmonic  coefficients 
of  the  geopotential  can  be  recovered  from  SST  data,  one  needs  a  model  that 
relates  the  information  in  thesedata  to  those  coefficients.  A  rigorous 
approach  involves  the  solution  of  many  variational  differential  equations 
for  the  two  satell ites. which  could  be  done,  in  an  average  sense,  by  the 
"analytical"  method  so  well  described  in  Kaula's  "Satellite  Geodesy"  (1966), 
and  in  an  instantaneous  sense,  by  the  "numerical"  approach,  an  example 
of  which  is  the  theory  behind  the  famous  "GEODYN"  computer  programme  (Martin 
et  al - ,  1970).  Regrettably,  except  for  some  major  break-through  in  computer 
science,  the  exact  application  of  either  technique  to  the  global  study 
of  SST  data  from  low-orbiting  satellites  does  not  seem  feasible.  The  problem 
is  the  sheer  size  of  the  spherical  harmonic  model  needed  to  represent  this 
data  realistically,  with  N>  300,  as  suggested  in  the  previous  paragraph, 
and  some  4xl06  observations  over  a  six  month  mission.  The  largest  models 
obtained  to  date  (for  instance  GEM  9,  by  Lerch  et  al.,  1977)  by  either 
technique  from  satellite  tracking  data  alone  have  not  exceeded  N=30,  and  have 
already  involved  very  lengthy  operations.  This  does  not  mean  that  these 
methods  have  no  role  in  the  analysis  of  SST  data:  the  "numerical",  at 
any  rate,  has  been  used  already  for  the  recovery  of  gravity  anomalies  from 
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Apollo-Soyuzs  data  (Kahn,  1979),  and  for  the  error  analysis  of  a  GRAVSAT 
mission  (Douglas  et  al.,  1980),  but  all  this  has  been  done  on  a  local  basis, 
and  the  purpose  of  this  report  is  to  look  at  the  problem  globally. 

If  a  rigorous  approach  is  not  practicable,  then  one  must  seek  some 
reasonable  approximation  that  makes  the  task  easier.  In  this  work  the 
author  has  followed  the  example  of  previous  studies  (Hajela,  1978;  Runnel 
1980),  assuming  that  the  time-derivative  of  the  line  of  sight  velocity 
can  be  approximated  to  a  sufficient  extent  by  the  line  of  sight  component 
of  the  inertial  acceleration.  As  explained  in  the  previous  paragraph, 
the  inertial  acceleration  due  to  gravitation  has  an  average  component  along 
that  line  that  is  not  zero,  due  to  the  powerful  zero  harmonic  GMr"2  and, 
to  a  much  lesser  extent,  to  the  even  zonals  (consider  (1.12, a)  with  m=0). 

A  non-zero  mean  acceleration  would  bring  the  two  satellites  (in  this  case) 
together,  contradicting  the  assumption  that  they  can  follow  (with  the  right 
initial  conditions)  the  same  orbit  with  the  same  mean  angular  velocity. 

The  effect  of  the  zero  harmonic  alone,  for  a  height  of  160  km  and  p=300  km,  is 

do  -  2  sin  (-^-)  s  43  gal . 

and  the  interpretation  of  this  is  plain  enough:  in  the  field  of  a  central 
point  mass  two  objects  initially  at  rest  would  fall  along  the  lines  joining 
their  initial  positions  to  the  attracting  mass  and,  because  their  motions 
converge  at  this  point,  they  would  be  moving  closer  to  each  other  as  well. 
Clearly,  this  is  not  the  case  when  the  two  bodies  turn  in  the  same  circular 
orbit,  where  the  relative  velocity  and  its  derivative  are  always  zero. 

The  approximation,  it  can  be  argued,  can  be  much  better  for  higher 
frequency  effects,  i.e.,  the  departures  of  both  velocity  and  acceleration 
from  their  values  for  a  central  mass  field  caused  by  the  mass  anomalies. 

The  question  is  too  complex  to  be  settled  by  a  simple  argument, 
so  the  conclusions  arrived  at  by  previous  authors  and  by  this  on^  too, 
for  that  matter,  are  of  necessity  no  more  than  educated  guesses  of  a  pro¬ 
visional  nature,  in  the  discussion  that  follows,  I  start  by  deriving  rigorously 
the  relationship  between  velocity  and  acceleration,  and  then  introduce 
"order  of  magnitude"  estimates  for  some  of  the  factors  involved. 

To  get  directly  to  the  results  of  interest,  one  should  consider  the 
residual  line  of  sight  velocity 

<5v  12  -  v*2  -  ijl  (1.20) 

where  is  the  velocity  computed  by  integrating  the 

equations  of  motion  numerically  with  the  reference  potential  model  U  : 
this  can  be  called  the  reference  velocity.  The  residual  velocity  is  what 
is  likely  to  constitute  the  data  in  a  real  life  situation,  and  has  been 
taken  for  such  in  the  error  analysis  of  Section  3.  Removing  the  reference 
velocity  (at  least  in  the  ideal  case  where  the  model  truly  represents  the 
field  up  to  degree  and  order  NM )  has  the  advantage  of  taking  away  the  very 
large  effect  of  the  second  zonal  whose  presence  is  undesirable  from  a  numerical 
point  of  view.  As  the  computed  orbit  is  likely  to  reflect  very  closely 
the  effect  of  the  attraction  of  the  Sun,  the  Moon,  and  of  the  major  planets, 
as  well  as  the  indirect  effect  of  the  tidal  deformation  of  the  Earth,  the 
residual  data  are  likely  to  be  freer  from  these  unwanted  gravitational  signals 
than  the  original  data,  resulting  in  less  interference  with  the  desired 
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information.  Non-gravitational  effects  are  removed,  presumably,  by  the 
drag-compensating  mechanism  . 


The  time-derivative  of  the  residual  line  of  sight  velocity  is 


Sviz  =  a!  (xLe^-xIi^iib  =j&2  ei2  -  £i2^Tei2^xI2  £12  -  xj^Te|^ 


=  xla  e12.x^)l(1c2)+il2(ix2-xI2ei2£x2)P^ *-kiC2)T{&)-k{i)%£W 


■l(c) 


=[xi2- xS^^eiz  +  p-1  [xl2  xi2-  x{‘>T +  (x<f>T e, -  (_xl2ei2)2] 

(1.21) 


assuming  that  the  computed  orbit  differs  from  the  true  orbit  by  only  a  few 
meters,  as  it  is  believed  to  be  the  case  when  using  contemporary  tracking 
data  and  field  models,  so  £i2sei2  and  |SC«  p  .  The  first  term  in  (1.21) 
is  the  residual  inertial  line  of  sight  acceleration,  so  the  second  term 
corresponds  to  the  discrepance  between  <5vi2  and  6ai2  =  ai2-a$4  •  This 
second  term  can  be  written 

e  =  P'1[(xl2  £n)2  -  (x^)T  en)2]  (1.22) 

where  en  is  the  unit  vector  pointing  outward  along  the  normal  to  the 
line  of~sight.  Therefore 

<5 v  1 2  ~  <$ai2  e 

*  6ai 2  +  o_1  6 ( x.T 2  £n  )2  (1-23) 


So  far  the  rigorous  analysis.  How  large  can  e  be?  Taking  the  mean 
value  of  e  over  all  rotations 

MU)*  -  M{6(xI2  V2}*  p-MI  S/  -  l  S**)  (1.24) 

m=l  m=l 

where  SnT  ^and  ^ ^re  the  average  power  at  m  cycles  per  revo¬ 
lution  in  x[\  en  and  xvC)T  en  ,  respectively.  These  are,  probably, 
of  the  same  order  of  magnitude~as  x|2  e12  and  x)f^  §12  .  As  the  model  is 
restricted  to  soacial  frequencies  of  NM  cycles  per  revolution  or  less, 
it  should  be  SfJ>  ( c)=o  if  m  >  NM  .  Replacing  %,h  with  the  Sm  in  (1-24), 
for  m>NM  the  summation  should  add  up  to 


m>  NM 

(e> 


=  „-i 


1 


m=NM+l 


(1-25) 


where  e  symbolizes  the  contribution  to  e  from  frequencies  above 
NM  .  Comparing  the  mean  value  for  the  mth  component  of  e  to  the  root 
mean  square  value  for  the  corresoonding  component  of  6a  J2  ,  or  Pi, 
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(1.26) 


wf 


.-i 


Pi  6  i 

-Em:  s3^i— 

in  "uj*  di1 


if  p  =  300 km  and  wsl.2  x  10*3  .  In  spite  of  the  presence  of  ,  this 
result  is,  in  fact,  dimensionless.  From  Table  1.1  follows  that  forNM=20, 
for  instance,  the  right  hand  side  of  the  expression  above  is  5.  x  10"'  , 

Rummel  (1981)  has  found,  by  a  similar  reasoning,  that  the  ratio  of 
M  f  5  x.1 2 }  to  the  total  rms  of  &alt,  (with  NM  =20)  may  be  ^ss  than  1.5  x  IQ"7, 
which  may  be  in  agreement  with  expression  (1.26):  the  relative  error  de¬ 
creases  with  frequency  because  of  the  m-2  factor  anc*  of  the  fas*,  decay 
of  with  m  ,  so  the  total  ratio  should  be  less  than  the  part^a1  '■atio 
at  m=20  .  Clearly,  there  are  several  assumptions  involved  in  these  results 
which  have  no  obvious  justification,  except  that  they  are  not  too  unreasonable. 
As  the  question  of  "reasonable"  or  "unreasonable"  is  a  tricky  one,  this 
author  would  conclude  that,  at  this  point,  the  omens  look  favorable  for 
the  model  of  the  line  of  sight  signal  proposed  here,  but  detailed  studies 
should  be  done  to  verify  this  question  further.  It  is  not  just  a  matter 
of  deciding  how  to  conduct  an  error  analysis  of  a  SST  mission,  but  also 
whether  it  may  be  possible  to  find  an  algorithm  for  reducing  the  data  from 
such  a  mission  which,  along  the  lines  explained  in  section  5,  could  be 
capable  of  resolving  potential  coefficients  to  a  much  higher  degree  and 
order  than  it  is  possible  at  present  with  existing  techniques.  One  way 
of  conducting  a  more  conclusive  study  of  this  matter  may  be  to  compute 
simulated  orbits  using  a  field  with  a  broad  spectrum,  such  as  that  of  the 

point  mass  model  used  by  Wagner  and  Colombo  (1979),  making  this  model  rotate 

like  the  Earth,  and  to  compute  also  the  orbits  corresponding  to  a  "model 
potential"  consisting  of  the  first  NM  harmonics  of  the  point  mass  field. 

From  both  sets  of  orbits  one  can  obtain  all  the  information  needed  to  calcu¬ 
late  the  two  terms  in  (1.23)  and  to  compare  them  with  each  other  in  various 

ways,  thus  throwing  a  much  brighter  liqht  on  this  whole  subject.  The  field 
of  the  point  mass  model  referred  to  above  has,  roughly,  the  same  Dower 
spectrum  as  that  of  the  Earth,  but  it  is  very  easy  to  compute  because  it  con¬ 
sists  only  of  200  mascons. 


Not  only  the  issue  of  the  accuracy  of  the  model  proposed  here  for 
the  line  of  sight  signal  can  be  clarified  by  numerical  studies,  but  perhaps 
better  models  may  emerge  that  share  with  the  present  one  all  the  practical 
advantages  and  that  are  also  truer  to  the  real  situation. 

The  model  chosen  for  v12  in  this  work  is  not  ai2  but,  as  explained 
in  the  previous  paragraph,  a12  =  ai2  -  a0  ,  where  a0  is  a  term  independent 
of  p  and  X  ,  about  which  more  will  be  said  in  section  2,  and 
equal  to  the  contribution  to  a12  of  the  zero  harmonic  and  of  the  zero 
frequency  terms  in  the  Fourier  expansions  of  the  even  zonals.  It  could 
be  argued  that  the  difference  between  this.and  the  model  used  by  previous 
workers  (who  like  Rummel,  have  chosen  v12  =  a12)  is  trivial,  because 
a„  must  be  almost  entirely  removed  when  the  reference  velocity  is  substracted 
from  the  data.  This  is  true  only  if  the  orbit  along  which  the  reference 
velocity  is  computed  coincides  with  the  true  orbit.  In  general,  this  is 
not  the  case,  as  orbital  errors  are  present  in  the  reference  orbit.  Because 
of  these  errors,  the  "removal"  of  the  zero  harmonic  and  even  zonal  effects 
accomplished  by  substracting  the  reference  velocity  is  not  as  thorough 
as  if  such  effects  had  never  been  included  in  the  definition  of  the  time- 
derivative  of  v i 2  ,  as  proposed  here.  If  left  in,  the  very  large  pueudo-effect 
of  the  zero  harmonic  would  give  the  impression  that  orbital  errors 

-15- 


2.  The  Mathematical  Model 


This  section  presents  the  observation  and  normal  equations  for  the  adjust 
ment  of  spherical  harmonic  coefficients  of  the  Earth's  gravitation  potential 
from  low-low  tracking  data.  These  equations  are  derived  here  under  some  sim¬ 
plifying  assumptions.  The  admissibility  of  such  simplifications  is  discussed 
elsewhere,  particularly  in  sections  1  and  4. 


2.1  Simplifying  Assumptions 

Computations  involving  spherical  harmonic  expansions  can  be  speeded 
up  greatly  if  there  are  regularities  in  the  distribution  of  the  data.  While 
such  regularities  may  not  occur  in  reality,  actual  and  ideal  distributions 
may  be  close  enough  to  each  other  to  ensure  that  the  results  obtained  for 
the  ideal  situation  are  also  applicable  to  the  real  one.  The  assumptions 
in  question  are: 

1)  both  satellites  describe  coplanar,  circular,  polar  orbits  with 
the  same  geocentric  radius  R  and  with  the  same  mean  angular  velocity  w  ; 

2)  the  plane  of  the  orbits  is  fixed  in  inertial  space,  and  fluctuations 
in  the  geocentric  angle  ^  between  the  two  satellites  are  disregarded; 

3)  the  Earth's  angular  velocity  vector  is  fixed  in  inertial  space, 
its  magnitude  is  constant  and  equal  to  ft  ,  and  its  direction  coincides 
with  that  of  the  figure  axis; 

4)  the  mission  lasts  an  integer  number  of  days  Ng  ;  oj  and  ft 

are  commensurable,  so  the  groundtracks  of  the  satellites  repeat  themselves 
(for  the  first  time)  after  Njg  days;  the  total  number  of  revolutions  of 
the  mid-point  between  satellites,  Nf  ,  is  prime  with  resDect  to  Nq  ; 

5)  there  is  perfect  compensation  for  ionospheric  propagation,  for 
radar  pointing  errors,  and  for  all  non-gravitational  forces  such  as  drag 
and  solar  pressure;  the  attractions  of  the  Sun,  Moon  and  major  planets 

have  been  accounted  for  exactly  when  computing  the  satellite  reference  orbits 
and  velocities  (the  data  consists  of  residual  velocities,  i.e.,  differences 
between  measured  and  reference  values); 

6)  data  are  sampled  at  constant  intervals  of  At  seconds  without 
interruptions  during  the  whole  mission;  there  is  an  exact  number  Np 

of  sampling  intervals  in  the  total  time  T  =  Nq  x  24  x  3600  s  ,  and 
Np  =  is  an  even  number;  ^  .  ..  ,  ... 

7 fz  the  sampled  values  consist  of  the  residual  line  of  sight  velocities 

averaged  over  Aa  seconds,  where  Aa  <■  At  ; 

8)  all  data  errors  are  uncorrelated,  have  zero  mean,  and  the  same 

standard  deviation;  , 

9)  the  line  of  sight  inertial  acceleration  differs  from  the  time  deriva¬ 
tive  of  the  line  of  sight  velocity  only  by  a  function  of  r  (constant 

for  a  circular  orbit);  • 

10)  at  satellite  altitude,  the  detectable  gravitational  signal  is  confined 
to  harmonic  terms  of  degrie  no  larger  than  N  =  331  ;  the  highest  frequency 

in  such  terms  is  less  than  half  of  the  sampling  frequency  ^  c/s  ; 
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11)  The  only  source  of  uncertainty  is  the  presence  of  errors  in  the 
measured  line  of  sight  velocities. 


Assumptions  (5).  (6),  (7),  and  (8)  by  and  large  state  what  a  perfectly 
successful  mission  should  produce  in  terms  of  data;  assumptions  (9)  and 
(10)  have  been  explained  already  in  section  1;  assumption  (11)  is  basically 
sound  if  assumption  (9)  is  admissible,  because  as  mentioned  in  paragraph 
(1.3)  and  further  argued  in  section  4,  (9)  implies  that  the  coupling  between 
orbit  determination  and  field  model! ing  is  weak,  so  orbital  errors* 
provided  they  do  not  exceed  a  few  meters,  have  little  effect  on  the  recovered 
potential  coefficients.  Assumptions  (1)  through'  (4)  define  a  simplified 
geometry,  the  validity  of  which  is  treated  in  detail  in  section  4. 


2.2  The  Extended  Legendre  Function 
The  equation 

m 


"nm 


U)  .  (?n!)  cos”*  [Slnn-J  -  1  sin"-”}2 

*  2nn! (n-m) !  2(2n-l) 


*  (n-»i....lnrm-3)  si„n-m-4  _  _  _  j 
2.4(2n-l)(2n-3) 


(2.1) 


defines  the  function  Lnm  of 

0 

that  has  the  following 

properties: 

(a) 

l„m(o)  =  p„m(siniJ>) 
nm*v/  nm*  ' 

and 

dh 

lmM  =  Pnm(sino) 

if  t 

(b) 

Lnm^  ”  Lnm^_(J>- 

if 

n  - 

m  is  even 

(2. 2, a) 

*"nm^  ~*'nm^’<^ 

if 

n  - 

m  i s  odd 

(2.2,b) 

Ln_  (o )  =  L  (fl'"4>) 
nm*  '  nm 

if 

m 

is  even 

(2.2,c) 

lnm(4>)  =-Lnm(^) 
nm*v/  nm 

if 

m 

is  odd 

(2.2 ,d) 

in  the 

interval  0  s  < 

2rr  . 

From  (a)  and  (b)  one  can  infer  that 

(1)  is  even  if  n-m  is  even,  odd  if  n-m  is  odd; 

(2)  L^U)  has  a  finite  Fourier  expansion  where  the  highest  term  is  the  nth 
harmonic  and  only  sines  or  cosines  are  present,  depending  on  the  parity  of 
(n-m): 


Lnm^  ^ 

I  C‘  cosp$ 
p=0 

if  n-m 

is  even 

(2. 3, a) 

Lnm(<f)  " 

r  ^msinP<p 

if  n-m 

is  odd; 

(2.3, b) 

P=0 


(3)  Lnm(d>)  has  half  wave  symmetry  (Lnm(p)»Lnm(TT-<f>))  if  m  is  even,  and 
half  wave  antisymmetry  (Lnm(l>)=-Lnm(1T"<i)))  if  m  is  odd.  Sums  of  sines 
or  of  cosines  which  have  such  symmetries  can  contain  only  even  or  odd  harmonics, 
so  the  Fourier  expansion  of  Lnn,  must  have  only  even  terms  if  n  is  even, 
and  only  odd  terms  if  n  is  ood.  As  a  consequence,  the  Fourier  coefficients 


hnm  with  p  of  opposite  parity  from  n  are  all  zero; 

(4)  tnm^)  "is  continuous  and  infinitely  differentiable  in  0  s  j  <  2- 
ana  takes,  together  with  all  its  derivatives,  the  same  values  as  Pnm  and 
all  its  derivatives  in  the  interval  T-<'h] t  • 

The  function  Ln„  is  the  analytical  continuation  of  Pnm  in  the  interval 
0  <  4>  <2tt  ,  and  can  be  called  for  this  reason  the  extended  Legendre  function  of 
the  first  kind,  degree  n  and  other  m  in  0  s  $  s  2n  .  Multiplying  Lnm 
by  the  same  normalizing  factor  as  Pnm  one  obtains  the  fully  normalized 
extended  Legendre  function  ~  ^ 


■ 


/2n+l 


LrOTi({,) 


l->> 


if  m  =  0 
otherwise 


(2.4) 


Consider  now  the  spherical  harmonic  of  degree  n  and  order  m 


(2.5) 


The  maximum  circle  containing  both  poles  and  the  point  (<j>=0,x)  on  the  equator 
consists  of  two  meridians,  of  longitudes  A  and  A+it  ,  respectively.  Along 
this  circle,  points  can  be  ordered  according  to  a  parameter  <j>’  in  the 
interval  0  <  <t>'  <  2rr  ,  which  increases  continuously  from  (q'=0,\)  towards 
the  N  pole  and  beyond.  The  geocentric  latitude  <j>  ,  on  the  other  hand, 
first  increases  towards  the  pole  like  <j>'  ,  but  on  crossing  the  pole  begins 
to  decrease  again  as  it  approaches  the  equator  at  (.j>=0,A+tt)  ,  and  it  is  negative 
in  the  southern  hemisphere.  The  harmonic  is  a  continuous  function 
of  $'  along  the  meridional  circle,  the  same  as  its  derivatives,  and  it 
follows  from  (2.5)  and  (2.2,c-d)  that 

C  ©  “<*«) 

■  (-1)"  Pnm  (sin*)  Q  n,  X 

*  Cn,n  <’"*>  ©  "  * 

■U*')  © 

This  relationship  shows  that  by  using  Cnm  and  q’  instead  of  Pnm  and  $ 
one  can  formulate  ?$m  avoiding  the  complications  that  would  arise  otherwise 
because  of  the  discontinuity  in  the  va’ue  of  the  longitude  at  the  poles. 

This  is  not  a  real  discontinuity  in  the  function  ,  but  merely  a  conse¬ 

quence  of  the  way  in  which  longitudes  are  defined. 

Replacing  '  with  its  Fourier  expansion  (2.3,a-b)  the  harmonic 
becomes 

Cn  =  I  cosre’  {J?*}  m  A  if  n  -  m  is  even 

p=0  p  s  (2.7) 

l  F£m  sin Pii*  m  A  if  n  -  m  is  odd 

p=0  w  s  1 " 
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r 


where  nm 

/Zn+T  h™  if  m  =  0 

^nm  _  P 

P  /2(2n+i)  in-m) T  h£m  otherwise 

fn+m) !  M 

and  there  are  orly  even  Fourier  terms  if  n  is  even,  and  odd  terms  if  n 
is  odd. 


2.3  Time  Series  Expression  of  the  Inertial  Line  of  Sight  Acceleration 

Reasoning  as  in  paragraph  (1.2),  but  considering  the  line  of  sight 
oriented  along  a  meridian  (polar  orbit)  instead  of  along  the  equator,  one 
gets  the  following  expression  for  the  line  of  sight  relative  inertial  accel¬ 
eration 

a12  *  a2  -  a2 

*  e[2  (rJSl)  ar  (S2)  +#15  a.  (SO  -  r(0Sz)  ar(S2)  +  aj$2) 

a  -(ar(Si)  +  ar (S 2 ) )  sin4j-  +  (a^Si)  -  a^(S2))  cos-j-  (2.8) 

where  terms  containing  the  "across  track"  acceleration  a,  have  been  dropped 
because  is  always  normal  to  the  orbital  plane  and  to  the  line  of  sight. 
Writing  the  radial  and  "along  track"  components  ar  and  a^  according 
to  (1.8,a-b),  replacing  Pnffl  with  Cnm  ,  and  choosing  the  mid-point  coordin¬ 
ates  <p'  =  +o-  ,  A  =  A i  +  A 2  as  independent  variables, 

2  ~2 
N  n 


ai  2 ( ^ » ’  »^) 


GM 

"R7  n=0  m=0 


*  W  cosr  i^r1 n™cos"J 

+  Snm  sinmA  ]  (2.9) 

where  N  is  the  smallest  degree  such  that  C°L  =0,  n>N  ,  according  to  the 
band-limited  assumption. 


Since 


L„„  <♦'  *  e)  ■  p^o  ™  ®  o(*'  *  B) 


(cosp(<j)"+S)  if  n-m  is  even,  sin p(cp+s)  if  n-m  is  odd),  it  follows  from 
elementary  trigonometric  relationships  that 

r  _  y  cnm  ,cos  pp'cos  P3  -  sin  p$'  sin  ps  . 

nm  ^  p^0  p  ‘sin  p<j>'cos  p@  +  cos  P4>‘  sin  pg  ! 


so 


lnm(*‘-1r)  +  Lnm^'+4)  "  2Jn  fT’c0S  p"£  ©  p ♦ 


n 

I 

p«0 
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A 


while 


Cnm(0''  '  Cnm^'+  =  2  n’nD  T  P*' 

p30 


and  therefore 


3T  [Cnm(0'-  T}  '  Lnm(4'+  T}  ]  =  2p£Q  Rp"  P  sin  P?"  S 


n  nm  _  .  i  i cos ^  (2 . 10 ,b) 

sm  p  J  .  _ 


From  expressions  (2.9)  and  (2.10,a-b)  follows  that 

GM  *  n  n 


ai2(R,  <P‘,  X)  =  [  I  C  I  2F™  (n+1)  cos  p  X  pi 1 

R  n=0  m=0  p=0  p  £  S1n 

+  ?  2fiJmosin  p.  cos-|  {cs^}  pd*]  (f)n  [Cnmcos  mA  +  $nm  sin  mX  ] 

p=0  (2.11) 

is  the  relationship  between  the  value  of  the  line  of  sight  relative  inertial 
acceleration  for  the  two  satellies_and  the  spherical  harmonic  coefficients  of 
the  gravitational  potential,  the  .  Carrying  out  various  multiplications 
indicated  in  (2.11),  and  making  use  of  the  trigonometric  equations 

2cos  Pi'  cos  mA  =  cos(pi‘+  mA)  +  cos(po’-  mX) 

2cos  d;1  sin  m\  =  sin(pi‘+  mX)  -  sin(pi’-  mX) 

2sin  pi'  cos  m\  =  sin(pi'+  mX)  +  sin(p$'-  mX) 

2sin  pi‘  sin  mX  =  -cos(pi'+  mX)  +  cos ( p<p ' -  mX) 


leads  to  tn.  expression 


N  n 


a.,(R  *•  X)  =  WJ_  y  V  rc  7  anm  |COS(p<i,+mX)+cos(p<!),-mX) 

ai2(R,  4>  ,  X)  ^  LCnm  (T)  ^  ap  lsin(p4).+mx)+sin(po'-m\) 

.  ?  /  r  nm  rsin(pi'+  mX)  -  sin(p4>'-  mX)  , 

+  nm  {W  4n  P  '-cos(p$'+  m\)  +  cos(p<t>’-  mX)* 


where 


a™  =  hpm  [(n+1)  cospf  sin|-+  P  sin  p|-  cos-|-] 


(2.13) 


so 


.nm 


if  n  and  p  have  different  parities  (par.  2.2,  remark  (3)). 


Introducing  time  as  the  independent  variable  through  the  formulas 
^  ”  ^ut^mod  2u 


and 


■[at] 


mod  2tt 


and  introducing  the  scaled  harmonic  coefficients 

r  =  p  t  a\n  GM 
Cnm  ^nm  '"TT'  r" 


nm 


*  ,a  ,n  GH 
Snm 


(2.14 , a) 
(2. 14,b) 

(2. 15, a) 
(2. 15,b) 
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expression  (2.13)  becomes 

,  _  r  r  ~  r  nm  jcos(pw  +  (itl)t  +  cos(pw-  rnR)t  , 

12(  '  '  n=0  m^Q  Cnm  p^0  ap  sin(pw  +  nfi)t  +  sin(pw  -  nfi)tf 

?  v  nm  r -sin (poo  +  nft)t  +  sin(poo  -  irtf)^  /- 

^nm  ap  'cos(pu)  +  irtl)t  -  cos(pu;  -  mfi)t  ‘  '  ‘  1 

This  last  formula  shows  au  as  a  time  series  representable  by  a  finite  Fourier 
series  where  the  angular  frequencies  present  have  all  possible  values  pw  ±  nil,  o< 
P  ,  m  i  n  ,  and  the  Fourier  coefficients  are  sums  of  terms  of  the  type  Cnm  a^m 
and  Snn  anm  ,  respectively .  The  time  origin  chosen  should  not  influence  the 
outcome  of  this  analysis.  The  choice  implied  by  (2.14,a-b)  corresponds  to 
t  =  0  when  the  mid-point  between  the  satellites  is  directly  over  the  equa¬ 
torial  point  (0,0). 


2.4  The  Correction  Term 


As  explained  in  paragraph  (1.2),  the  model  for  the  time  derivative  of  the 
relative  line  of  sight  velocity  used  in  this  study  is  not  the  relative  inertial 
line  of  sight  acceleration  ai2  ,  but  this  acceleration  minus  a  time-invariant 
term  a0  (expression  (1.13)) .  or  modified  acceleration  ai2  .  According 
to  (2.16),  this  time-invariant  part  can  only  be  due  to  terms  where 
Pu)  ±  mft  =  0  .  As  indicated  later  in  paragraph  (2.6),  one  of  the  assumptions 
made  in  this  study  imply  that  pw±mfi^0  if  0  <  m  £  N  so  this  leaves 
only  the  case  p=Q  ,  m=0  ,  corresponding  to  the  even  zonal s,  n=0  in  partic¬ 
ular.  Therefore,  the  corrected  acceleration  is 

a»:(t)  =  a12(t)  -  a<, 


where 


N 

■Jo 


no 


,  a  ,n  GM 

(TT)  -p 


(2.17) 


is  a  term  that  depends  only  on  R  ,  as  expected. 


2.5  The  Observation  Equation 

One  of  the  assumptions  made  in  paragraph  (2.1)  was  that  the  orbits 
are  periodic  with  a  period  T  equal  to  the  length  of  the  whole  mission. 
Consequently,  the  various  angular  frequencies  (pw±  mfl)  present  in  the  right 
hand  side  of  (2.16)  are  harmonics  of  fundamental  frequency  <^o  =  ^  .  The 
modified  line  of  sight  relative  inertial  acceleration  812  is  regarded 
here  as  the  true  time  derivative  of  the  relative  line  of  sight  velocity, 
in  accordance  with  expression  (1.13).  Therefore 

v 1 2 ( t )  =  J*  3i2  (t)  dt  +  vi2(0)  (2.18) 

0 


Because  of  the  assumption  of  orbital  periodicity,  replacing  the  integrand 
in  (2.18)  with  the  difference  between  the  right  hand  side  of  (2.16)  and 
a0  according  to  (2.17),  and  integrating  this  difference  term  by  term,  one 
gets  the  following  expression,  where  no  zero  frequency  term  is  present: 
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(2.19) 


Vl  2  (t ) 


I  I 


nm 


n 

IaP 
p=z  K 


nm 


sin(pu)  +  mfl)t  + 
pw  +  m  n 

■cos(pu  +  m$t  - 

pw  +  m  n 


sin(pu)  -  mft)t 
pu  -  m  J1 

cos ( pu  -  nfl)t 
pu)  -  m  n 


cos(p ui+m Sl)t  -  cos(pu-  mSl)t 

pu  +  m  u  p  oj  -  mil 

sin(pu  +  nfrt)t  -  sin(p^)  -  mA)t 
pu  +  mil  pu'  -  m* 


where 

z 


0  if  MO  (2.20) 

1  if  m  =  0  and  where  C13  ,  Cn  ,  Sn  have  been  dropped  as 
unknowns,  because  the  origin  of  coordinates 
coincides  with  the  geocenter. 


The  actual  observations  consist  of  time  averages  V12  of  the  measured 
relative  line  of  sight  velocity.  This  averages  are  taken  over  identical 
intervals  of  length  Aa  seconds,  spaced  At  seconds  apart.  Tne  signal  in 
the  observations  is,  then,  the  line  of  sight  velocity  averaged  over  Aa  : 


*”(t)  =  Ta  h-  Aav^(t)  dx 

,  N  n  n 

=  "Aa”  l  L  ^nm  J  “p 

n=2  m=C!  p=z 


nm 


S(  (pu+  nfllt,Aa)  +  S((pu  -  mft)t,A  a( 

(pu  +  mfl)2  (pu  -  mh)2 

-C((pu  +  mft)t,Aa)  C(  (pu  -  m£l)t,A  a) 
+  m$2 - - 


n 


+  Snm  l  ap 
P=z  H 


nm 


where 


CLLe^l^l^Aa)-  C(_(p_u-  ,.^)t,Aa) 

(pu  +  mfi) “  ( pu  -  mQ ) • 1 

S((pu  +  mH)t,Aa)-  S((pu-  mfi)t, Aa 
urrwp  (pu  -  m:)- 


(2.21) 


and 


C((Pu  t  nv2)t,  Aa)  =  cos(pu)  ±  mC)t  sin(pu  ±  nfl)  Aa 

+  sin(pu  ±  mn)t  (l-cos(pu  ±  mil)  Aa) 

S( (pu  ±  mfi)t,  Aa)=  sin(pu  ±  rtl)t  sin(pw  ±  mfi)  Aa 

-  cos(pu  ±  mfi)t  (l-cos(pu  ±  m£2)  Aa) 


(2. 22, a) 


(2.22, b) 


°^-s-erved-  Vd^Llg.  t.j  consists  of  v12  plus  the  measurement  noise 
average  over  tj-  aa  *  t  <  tj  ,  or  ni  :  - - 


v!5i}  +  M  =  vfo\obs) 


(2.23) 


Rearranging  the  order  of  summation  with  respect  to  n  and 
and  replacing  the  result  in  (2.23)  one  arrives  to  the  observation 
it  will  be  used  in  the  error  analysis  - — 


m  in  (2.19) 
equation  as 


N  N 

L  A  ^nm 

$ax(m,2) 


n 

l 


p=z 


.nm 

ap 


+  s 


n 

nm  £ 
p=z 


nm 


C((p,o  +  mn)tj  ,a  a)  -  C(U-mfl)ti,A  a) 
(p<ij  +  m£>)5  (Pa;  -  mo)2  1 

S((puj  +  mflltj  ,Aa)  -  S( (pu  -  mft)tj ,Aa) 1 
(pw  +  mft)2  [pLj  -  mfi)2  ) 


.( t,) 

"  V  (obs) 


ri 


(2.24) 


where  max(m,2)  indicates  the  largest  of  m  and  2  ,  and  where  r^  is 
the  residual  or  difference  between  the  measurement  and  the  value  calculated 


by  replacing  the  Cnm  with  numbers  in  the  left  hand  side.  When  these 
numbers  are  the  true  values  of  the  unknown  C^m  and  the  model  is  perfect 
(as  it  is  supposed  to  be  in  the  case  here),  then 


ri 


ni 


(all  this  is  in  keeping  with  established  practice  in  geodetic  literature) . 
Consider  the  following  vector  notation 


_  p- ( 0)  (At)  u  (*At) 

'  LVl2(obs)  Vl2(obs) *  •  •  Vl2  (obs)‘ 

(!-At«  (Np- 

- 1  )At)  -jT 

-1 2 (obs ) 

•  •  x"(obs) 

(2. 25, a) 

r 

=  [r0  rx  •  •  •  ^  .  .  .  rNp-1]T 

( 2 . 25 ,b ) 

fr-mm  ^m+l)m‘  ’  ’  CNm^ 

( 2 . 25 ,c ) 

=  ^mm^m+l)m"  ‘  ’  W 

(2.25,d) 

c 

=  CcJ  cj  sj  cj  sj  .  .  .  c^  s^  .  .  . 

cT  ST-,T 

( 2 . 25 ,e ) 

where  _Vi;[0bs)  and  -  are  bottl  NP  “  vectors  and  c  is  a  Nc  -  vector, 

Np  being  the  total  number  of  observations  (some  3.9  x  106  measurements 
over  six  months  if  At  =  4)  and  Nc  =  (N+lJ2-  3,  the  number  of  coefficients 
in  the  band(*)  2  <  n  s  N  .  The  set  of  all  observation  equations,  or  system 
of  observation  equations,  in  matrix  notation,  is 

A£=  v12(ob$)  +  r  (2.26) 

where  A  is  the  Np  x  Nc  matrix  of  the  observation  equations^  The  unknowns, 
according  to  (2.24),  are  the  scaled  potential  coefficients  Cftm  ,  instead 
of  the  actual  coefficients  Cftm  ,  which  are  the  ones  desired.  Once  the 
Chm  are  known,  however,  the  can  be  obtained  by  a  trivial  operation 

based  on  (2.15,a-b).  The  same  is  true  of  the  standard  deviations,  which 
are  the  quantities  relevant  to  this  error  analysis. 

The  observation  equation  (2.24)  presents  the  relationship  between  data 
and  unknowns  in  time.  The  corresponding  relationship  in  space  can  be  obtained 
directly  from  (2.24)  by  replacing  <f>*  and  X  for  t  in  accordance  to 
(2.14,a-b).  The  temporal  representation,  however,  makes  the  overall  treatment 
of  what  is  to  follow  simpler,  and  is  adopted  here  for  this  reason. 


'  ^Coefficients  with  n  >  NM  are  actual  field  coefficients;  with  n  s  NM 
(NM  being  the  highest  degree  in  the  reference  field  model  mentioned  in  para¬ 
graph  (1.2))  the  C$m  correspond  to  residual  coefficients:  the  differences 
between  the  actual  coefficients  and  the  coefficients  of  the  model,  if  the 
data  consists  of  residual  velocities,  as  assumed  in  section  3. 
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If  Earth  rotation  and  satellites  revolution  are  congruent  over  the 
length  of  the  mission,  the  total  number  of  revolutions  Nr  and  the  total 
number  of  days  Nq  are  both  positive  integers.  According  to  condition 
(4)  in  paragraph  2.1  they  are  also  relative  primes,  i.e.,  without  common 
factors  other  than  the  unity.  As  shown  here  this  property  rules  out  certain 
relationships  between  the  frequencies  pw  ±  mfi  in  (2.24),  the  absence  of 
which  simplifies  the  mathematical  treatment  of  the  error  analysis,  as  it 
will  be  explained  in  paragraph  2.7.  Taking  m  in  the  interval  -N  <  m  <  N  , 
pw  ±  mfi  can  be  written  more  simply  as  pm  +  mf>  .  The  'elationships  of 
interest  have  the  form 


Pui  +  mft  »  ju>  +  qn  (2.27) 

where  either  p  f  j  ,  m  /  q  ,  or  both.  The  question  as  to  whether  such 
relationships  are  possible  can  be  answered  in  two  parts: 


Case  where 


If  (2.27)  is  possible,  then 


m  -  q 
J  -  P 


_  _  m  u)p *  _  Nr 

51  fi  uio 1  Nn 


(2.28) 


where  oo0  is  the  fundamental  angular  frequency  of  the  mission:  = 

As  both  Nr  and  Nq  are  positive  it  follows  that  =  -fe-  so 


I "i  “  <ll  Nr 


|j  -  Pi  N, 


As  |m-q|  Nq  and  |j-p|  Nr  are  both  positive  integers  and  Nr  and  Nq 
have  no  common  factors,  it  must  be  | m -  q [  =  Nr  K  for  some  K  z  1  ,  so 


|m  -  q|  >  N, 


(2. 29, a) 


As  -N  s  m  s  N  and  -N  s;  q  <;  N  ,  it  follows  that,  under  the  band-limited 
assumption. 


|m  -  q|  <  2N 

As  a  consequence  of  (2.29,a-b) 


Nf  <  2N 


(2.29,b) 


(2.29,c) 


if  (2.28)  is  true.  In  the  present  case,  the  number  of  revolutions  over 
six  months  with  the  satell  ites  at  a  height  oT  160  km  exceeds  2900,  while  the 
maximum  degree  n  of  the  harmonics  in  the  band  considered  here  is,  according 
to  paragraph  (1.2),  N  =  331  .  This  means  that  2N  <  Nr  ,  contradicting 
(2.9,c)  and  thus,  (2.28).  Therefore,  for  the  given  satellite  height  and 
field  bandwith,  (2.27)  is  impossible  if  in  f  q  . 

(b)  m  *  q: 

In  this  case  (2.28)  is  never  true,  because  the  first  member  is  always 
zero  except  in  the  trivial  case  p  *  j  ,  where  it  is  indeterminate.  Therefore 


(2.27)  is  also  impossible  when  m  =  q  and  ,  thus,  whenever  either  p  ^  j 
or  m  f  q  apply. 


An  immediate  consequence  of  the  impossibility  of  (2.27)  under  the  assump¬ 
tions  of  paragraph  2.1  is  that 

pw  ±  mft  =  0  (2.30) 

is  also  impossible,  except  in  the  trivial  case  p  =  0  ,  m  =  0  ,  which  has 
been  excluded  by  the  way  in  which  the  modified  line  of  sight  acceleration 
is  defined  (see  expressions  (1.13)  and  (2.17)).  Since  the  rate  of  precession 
of  the  node  of  a  polar  orbit  is  nil,  condition  (2.30)  corresponds  to  what 
is  known  in  satellite  geodesy  as  a  resonant  orbit.  So  the  assumption  that 
Nr  and  Np  are  relative  primes  excludes  resonances.  This  means  that  the 
orbit  is  notallowedto  repeat  itself  after  a  number  of  revolutions  that 
is  a  submultiple  of  Nr  ,  so  no  sub-cycles  occur  within  the  grand  cycle 
whose  period  is  the  length  of  the  whole  mission,  i.e.,  T  seconds  or  No 
days.  Because  of  this,  the  sampling  of  the  gravity  field  over  the  face 
of  the  Earth  is  the  most  even  that  is  possible,  as  with  resonances  many 
arcs  would  be  superimposed,  so  the  "footprint"  of  the  mid-point  between 
the  satellites  would  cover  the  ground  rather  coarsely.  Without  resonances, 
the  arcs  are  all  different  and,  thus,  better  spread  out. 

In  reality,  resonance  is  most  unlikely  to  occur  exactly:  the  actual 
situation  is  much  too  complicated.  How  realistic  is,  even  so,  the  assumption 
that  Nr  and  Np  are  relative  primes?  While  perfect  congruence  (both 
numbers  integers)  is  also  quite  unlikely,  the  situation  need  not  be  too 
different  from  that  implied  by  the  assumption.  Choosing  Np  =  179  (a  prime 
number)  and  Nr  =  2933  (another  prime),  both  are  then  relative  primes. 

For  the  satellites  to  orbit,  at  160  km,  2933  times  around  the  Earth  in  precisely 
179  days,  a  change  in  w  of  only  -4  parts  per  thousand  from  its  actual 
value  is  needed  or,  equivalently,  an  increase  in  ft  of  4  p.p.  thousand. 


2.7  The  Scalar  Product  of  Two  Columns  of  the  A  Matrix 


If  the 

of  observat,ionrequations 


Np  -  vectors 


inm 


and  Cj<q,  respectively,  then  their  scalar  product 

oai3  =  (aa  )T  a3  =  Y1  aa(ti)  a3'*^ 
pnmkq  '^nnr  ^kq  T  anm  akq 

i=0 


and  a£q  are  the  columns  in  the  matrix 
of  (2.26J  corresponding  to  the  scaled  coefficients 


(2.31) 


should  depend,  according  to  (2.24),  on  what  the  mission  parameters  Aa  ,  At  , 
and  i p  are,  and  also  on  the  values  of  sumsof  products  of  cosines  and  sines 
of  the  various  arguments  (pw  ±  mft)t-j  ,  where  the  discrete  values  t^  cover 
the  whole  length  of  the  mission.  The  products  of  the  columns  of  A  have 
to  be  obtained  as  part  of  the  formation  of  the  normal  matrix  of  the  adjustment, 
the  inverse  of  which  provides  the  a  posteriori  accuracies  that  are  the  main 
objective  of  the  error  analysis.  Sums  of  products  of  sines  and  cosines 
sampled  at  regular  intervals  are  strongly  influenced  by  the  relationship 
between  the  sampling  frequency  (inverse  of  the  sampling  interval  At)  and 
the  highest  frequency  in  the  arguments  of  the  trigonometric  functions  involved. 
In  particular,  it  is  important  to  know  whether  the  highest  frequency  is 
below  the  Nyquist  frequency  (one  half  the  sampling  frequency)  or  not,  in 
order  to  choose  the  most  convenient  treatment  for  those  sums  of  trigonometric 
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products.  In  the  case  at  hand  At  =  4  s,  so  the  Nyquist  frequency 
is 

"y  -  73J  cycles/s 
=  0.125  c/s 

while  the  highest  frequency,  corresponding  to  the  term  (Noo  +  NO)t  is, 
for  N  =  331. 


fmax  =  TH>  +n)  =  0-066  c/s 


Therefore, 

Wy  <2-32> 

or  the  highest  frequency  is  less  than  the  Nyquist  frequency.  Under  this 
condition,  the  following  formulas  apply: 


l  cos(pa)  ±  m^)ti  cosCju  ±  qfi)ti  =  £  sin(pu>  ±  mfljti  sin(ju>  ±  qSl)t, 
i=0  i=0 


||p-  if  (po  ±  ft)  ■  (jw  ±  qft) 


(0  otherwise 

p=j=m=q=0  being  excluded,  and  Np  being  even; 


(2. 33, a) 


V  cos(pu)  ±  msi)t.  sin(joj  ±  qfl)t.  =  0  always 
1=0  1  1 


(2.33,b) 


From  these  formulas,  and  from  expressions  (2.22,a-b),  one  gets 


l  C  ((pto  ±  mft)ti,Aa)  C  ((ju>  ±  qftJti.Aa) 

i=0 

,  ,  ,,  N  (1-cos  (pu>±mft)Aa)  if  pw±mn 

=  .1  S  ((pu)±mft)ti  ,Aa)S((jw+qft)t-j,Aa)  =  P  =  jui+qft 

1=0  0  otherwise  . 

(2. 34, a) 


p  =  j  =m=q*0  being  excluded,  and  Np  being  even; 

Np-1 

T  C  ((pw  ±  mnJtJ.Aa)  S  ((jw  ±  qft)t<),Aa)  «  0  always 
i*0  1 


(2.34,b) 


As  a  result  of  expressions  (2.34,a-b)  above,  of  the  fact  that  apm  *  0  if 
p  has  different  parity  from  n  ,  and  of  the  fact  that  pw  ±  mft  ?  jw  ±  qn 
if  at  least  p  f  j  ,  as  explained  in  the  previous  paragraph,  replacing 
and  aj^lti)  in  (2.31)  by  their  respective  expressions  according 
to  (2.24),  i.e.. 


aa(ti)  =  -L  ?  anm 
nm  p=z  P 


i  1 


lOi 


where  {  l"1  is  the  first  or  the  second  bracketin  (2.24),  depending  on  ot  , 
and  similarly  for  ajH^i)  ,  one  gets 


.even 


nmkq 


}  and  k  |oc^  j 

1  ‘even1 


h 

Aa 


if  a  t  6  ,  m  f  q  ,  n  {0(Jd 
irnn(n,k)^nm  cos(pjj-nnfl)  Aa)  +  (1-  cos(pa)-mil)  Aa)  ^ 

p  p  (pu+  mQ)1*  (pw  -  mft)4 


P=2 


min(n,k)  being  the  smallest  of  n  and  k  . 


(2.35) 


as  the  expression  for  the  scalar  product,  which  shows  that  under  the  assump¬ 
tions  the  product  is  zero  in  the  majority  of  cases.  This  results  in 
many  elements  in  the  normal  matrix  being  zero  as  well,  which  is  a  major 
advantage  when  setting  up  this  matrix,  as  only  the  relatively  few  non¬ 
zero  elements  have  to  be  computed.  Moreover,  as  shown  in  paragraph  (2.9), 
a  suitable  ordering  of  the  unknowns  groups  this  non-zero  elements  in 
a  block-diagonal  structure,  so  inverting  the  very  large  normal  matrix 
reduces  itself  to  inverting  the  much  smaller  diagonal  blocks,  a  crucial 
fact  as  far  as  the  feasibility  of  obtaining  the  a  posteriori  covariance 
matrix  is  concerned.  Finally,  expression  (2.35)  shows  that  the  non¬ 
zero  elements  can  be  computed  from  the  values  of  the  aGm  ,  which  are 
the  Fourier  coefficients  of  the  Lnm  multiplied  by  scale  factors  (expression 
(2.13)).  Because  the  apm  are  zero  when  p  and  n  have  different 
parities,  because  Lnm  are  expansions  of  sines  only  or  of  cosines  only, 
and  because  n  cannot  be  larger  than  N  =  331,  the  total  number  of  terms 
in  the  summation  of  (2.35)  is  much  less  than  that  of  the  terms  in  the 
summation  of  (2.31),  resulting  in  considerable  economies  when  computing 
the  non-zero  elements  of  the  normal  matrix. 

The  Fourier  coefficients  can  be  obtained  by  computing  [—,{<!>')  at 
regular  intervals  A<J>'  =  ■&-  (the  highest  frequency  in  Lnm  (V)  is 
that  of  cos  ncp '  or  sin  n$'  )  and  then  carrying  out  a  numerical  Fourier 
analysis,  or  discrete  Fourier  transform,  of  these  values.  This  is  done 
in  the  program  of  appendix  B  by  means  of  a  mixed-radix  Fast  Fourier  Transform) 
algorithm  chosen  on  grounds  of  efficiency.  Finding  the  required  value  of  Lnm 
in  the  interval  0  £  <p *  <.  2tt  is  simplified  by  the  use  of  expressions  (2.2,a-c) 
and  of  the  relationship  Cnm(4>')  =  Pnm(4>')  if  0  s  <J>‘  s  ■§■>  ,  which  reduces 
most  of  the  effort  to  that  of  computing  the  values  of  Pnm  at  regular 
intervals  in  0  <  o'  <;  . 


2.8  Least  Squares  Adjustment 
Consider  the  quadratic  form 

Q  =  rT  Pr  (2.36) 

where  P  is  a  syircnetric  N_  x  Np  matrix  and  where  r  is  the  vector 
of  residuals  p  H 


I  =  Ac  -  v12(obs) 
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(2.37) 


according  to  (2.24),  (2.25,b)  and  (2.26).  The  form  (2.36)  is  a  function 
of  £  ,  the  vector  of  potential  coefficients,  through  (2.37).  The  vector 
£  that  minimizes  the  quadratic  form  must  satisfv  the  condition 


<-§§-)T  •  (ATPA)  c  -  ATP5>Mobs)  =  0 

(where  0  is  a  null  Nc  -  vector)  so 
c  =  (ATPA)“  ATP^(0bs) 

and  the  form  has  a  minimum  at  c  provided  P  is  a  positive  matrix. 
Expression  (2.39)  above  can  be  written 

-  p  --  vobs) 


(2.38) 


(2.39) 


(2.40) 


where  Fp  is  the  optimal  estimator  matrix  corresponding  to  the  weights 
matrix  P  .  In  general  £  f  c  ,  because,  even  if  the  signal  in  the  data 
(v12)  and  the  unknown  parameters  (c)  are  related  exactly  to  each  other 
by  the  matrix  equations  vi2  =  Ac  ,  the  data  _Vi2(0bs)  contains  noise 
n  in  addition  to  the  signal.  The  noise  £  propagates  into  the  estimate 

£  =  Fp  ii2(obs)  =  Fp  (^12  +  ^  =  Fp  iiz  +  Fp  H  (2.41) 

resulting  in  a  difference  between  £  and  c  ,  or  error, 

£  -  Fp  is-  Fp  n  =-e^  (2.42) 


The  variance-covariance  matrix' of  these  errors  in  c  ,  or  a  posteriori 
random  errors,  is 

E  =  E  {e  e  T} 
n  _  i 


(2.43) 


where  E  {  |  is  the  usual  mathematical  expectation  operator  of  statistics. 

The  diagonal  elements  of  the  error  matrix  are  the  a  posteriori  variances  of 
the  estimated  coefficients,  or  formal  accuracies  of  the  adjustment.  Tf 

P”1  =  D  =  E  {n  nT}  (2.44) 


where  D  is  the  Np  x  Np  symmetrical  and  positive  matrix  corresponding 
to  the  data  errors  and  known  as  the  a  priori  variance-covariance  matrix, 
then 

En  *  MfD  "  (F0  £)T1  =  E  {FD  n  nT  F^ } 

*  Fo  Mu  isJ }  fdt  *  fd  0  fdT 

*  (At  D_1A)_1AT  d_1d  d_1a  (AT  O^A)-1 

-  (At  D^A)'1*  G'1  (2.45) 

where  Fn  *  (A^  0’IA)'1A^  0*1 
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Matrix 


(2.46) 


G  =  AT  D_1A 

is  known  as  the  normal  matrix,  because  (2.38)  can  be  written 

Gc  =  b  (2. 47, a) 

with 

b  =  AT  D'1  v12(ob$)  (2.47 ,b) 

and  (2. 47, a)  is  the  system  of  the  normal  equations.  Therefore,  when 
P  =  D_1 ,  the  a  posteriori  variance-covariance  matrix  is  identical  to 
the  inverse  of  the  normal  matrix.  In  particular,  when  all  data  errors 
are  uncorrelated  (E  (n-j  nj}  =  0)  and  all  have  the  same  standard  deviation 
c  ,  then 

D  =  o2  I 


where  I  is  the  unit  Np  x  Np  matrix,  so 
G  =  AT  a2  I  A  =  d2  At  A 


(2.48) 


and  the  elements  of  the  normal  matrix  have  the  form 


ga^, 

anmkq 


=  o 


•  :/aa  \T  / 
'^finr  -kq 


.  -2 
■  u 


a6 

pnmkq 


(2.49) 


where,  in  the  case  under  study,  p“£.  can  be  calculated  from  (2.35). 

Among  the  important  properties  of  the4least  squares  estimate  c  is  that 
of  being  a  minimum  variance  estimate  when  P  =  D"1,  which  means  that 
the  diagonal  elements  of  G-1  (and,  then,  their  sum  or  trace  tr  { En } ) 

are  minimized.  If  the  probability  distribution  of  the  error  is  gaussian, 
then  the  estimate  c  is  best  in  the  sense  that  the  a  posteriori  variances 
are  the  smallest  for  all  estimates,  linear  or  nonlinear.  Moreover,  the 
most  likely  value  of  c_  coincides  with  c  ,  so  £  is  also  the  maximum 
1 ikel ihood  estimate.  When  the  model  Vi2  =  Ac  is  perfect,  as  assumed 
here,  and  E  {nj  =  0  ,  then 

E  {£}  =  (AT  D"1Ar1A’  D_1(E  {v12}+  E  {n})  =  (AT  D^AfA1  D_1A  c  (2  g 

For  this  reason,  the  estimate  c  is  called  unbiased.  Expression  (2.50) 
further  indicates  that,  in  the  absence  of  noise,  the  estimates  are  identical 
to  the  true  values  of  the  unknowns. 


The  rather  impressive  list  of  properties  of  the  least  squares  estimator 
(2.39)  with  P  =  D-1,  together  with  the  relative  simplicity  of  the  theory 
and  rather  straightforward  nature  of  the  calculations  involved,  have 
made  this  type  of  estimator  the  most  widely  used  in  geodesy  as  well  as 
in  many  other  branches  of  technology  and  of  natural  science.  In  reality, 
models  are  never  exact,  statistics  never  truly  gaussian  (this  would  include 
cases  where  the  error  is  extremely  large,  but  such  occurrences  are  usually 
edited-out  during  pre-processing),  or  even  truly  known,  so  the  practical 
implications  of  the  theoretical  properties  of  the  estimator  are  somewhat 
obscure.  What  is  clear  is  that  (2.39)  with  P  *  D_1  is  a  "sensible"  way 
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of  processing  data,  as  it  gives  larger  weights  to  the  better  measurements 
(those  with  smaller  variances),  and  less  weight  to  the  worst  (larger 
variances),  at  least  when  0  is  diagonal.  If,  as  in  the  present  case, 
all  measurements  are  assumed  to  be  equally  good  (same  <?) ,  then  they  are 
all  weighted  equally.  From  the  point  of  view  of  the  error  analysis  which 
is  the  purpose  of  this  stuay,  the  determination  of  the  formal  accuracies 
requires  creatinaand  inverting  the  normal  matrix  G  . 


2.9  The  Structure  of  the  Normal  Matrix 


One  of  the  assumptions  made  in  paragraph  (2.1)  was  that  the  errors 
in  the  data,  the  elements  n-j  of  £  ,  were  random,  uncorrelated,  and 
had  all  the  same  standard  deviation  c  .  Consequently,  the  elements  of 
the  normal  matrix  G  can  be  calculated  from  (2.49)  and  (2.35)  which, 
combined,  give 


qa6 

-nmkq 


0  if  m  f  q  ,  a  /  8, 
Min(n,k) 


"CfNnd  k 


lodd 


rOdd 

‘even- 


,.2  N 


Aa- 


l 

P=Z 


nm  ,km 
P  P 


r(  1-cos (pm  +  mC)Aa)  (l-cos(p<d-  mfl)Aa-i 
L  (pu  +  mO)1*  (Du  -  m PI u  -1 


(2.51) 


As  already  pointed  out  (par.  (2.7)),  there  are  many  elements  in  G 
that  are  zero  and  it  is  possible  to  arrange  the  unknown  Cqm  so  that 
G  exhibits  a  very  convenient  structure.  Consider  the  ordering  given 
in  expressions  (2.25,c-e)  to  the  elements  of  c  .  where 


cT  sT 


C  =  [clcj  sjcjsl  . 

^m  ”  ^mm  ^(n+i)m^+2)m-  *  * 


•  4  %T 


-^m  ^mm  %i+l)m 


S  1 
‘  ’  •  inl¬ 


and  now  separate  each  Cj„  and  %  in  two  halves  cl  °,  cjrf1,  a"d  s , 


,  so  c^=0,  Stff0  contain  onTy  C'nm  where  n  is  even. 


-Y=* 


y=o 

-2®_,  >  , 

sl~l  only 


-Sj^i  9  ^JTl  9  _£jn  v-unwam  uuijr  mici  c  n  to  circn,  Ujp  j 

Cftm  where  n  is  odd.  With  the  unknowns  arranged  in  this  way,  the  cpi  , 
sm  contain  coefficients  of  the  .same  order  m  ,  the  Cnm 
"cosine"  coefficients  and  the  Srim  "sine"  coefficients,  and  each  of  these 
are  split  further  according  to  parity.  If  all  this  is  done,  then  G  can 
be  partitioned,  as  shown  in  figure  2.1,  into  blocks  Gmq  .  Those 
along  the  main  diagonal,  or  G™,  ,  correspond  to  unknowns  of  equal  order 
m  ,  and  are  further  partitioned  each  into  four  blocks,  according  to  a  and 
8,  and  each  of  these  four  blocks  into  another  four,  according  to  the  parities 
of  n  and  k  .  All  Gmn  ,  except  for  those  along  the  main  diagonal. 


In  the  diagonal  G^  only  the 


are  zero  matrices  according  to  (2.51). 

sub-blocks  where  a  =  3  and  n  and  k  have  the  same  parity  (dashed 
in  the  picture)  can  contain  nonzero  elements.  This  sub-blocks  are  diagonal  sub- 
blocks,  so  the  whole  normal  matrix  is  zero  except  for  those  diagonal  sub- 
blocks:  G  is  block  diagonal . 

More  concretely,  G  contains  4x(N+l)-4  diagonal  blocks  that  are 
not  null  matrices.  According  to  (2.51),  the  elements  of  "cosine"  blocks 
(a  *  8=  0)  are  identical  to  those  of  "sine"  blocks  (a  =  8=  1)  in  corres¬ 
ponding  positions.  This  means  that  only  the  2x(N+l)-l  "odd  n"  and  "even  n" 
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Figure  2.1:  structure  of  the  normal  matrix 


blocks  have  to  be  formed  when  setting  up  G  ,  and  only  they  have  to  be 
inverted  when  inverting  G  .  The  smallest  blocks  to  be  inverted  have  only 
one  element  and  correspond  to  m  =  N  ;  the  two  largest  blocks  have 
i(N+l):  elements  each  and  correspond  to  m  =  0  (zonals)  with  n  even 
and  odd,  respectively.  Altogether,  the  2x(N+l)-l  blocks  contain  in 
the  order  of  ^(N+l)3  different  elements,  as  the  blocks  are  also  symmetrical. 
In  the  case  of' a  general  symmetrical  matrix  the  total  number  of  different 
elements  to  be  computed  is,  approximately,  i (dimension)2  ,  or  i(N+l)“ 
in  the  case  at  hand.  The  reduction  in  calculations  when  setting  up  G 
is,  therefore,  of  the  order  of  ^-(N+l)  .  Furthermore,  expression  (2.13) 
shows  that  the  a^m  can  be  computed  from  the  scaled  Fourier  coefficients 
hpm  of  the  Lnm  .  As  many  hp*11  are  zero,  as  pointed  out  in  paragraph 
(2.7),  the  number  of  operations  using  (2.51)  rather  than  the  general 
expression  (2.49)  is  reduced  by  at  least  an  order  of  magnitude.  All 
these  savings  in  computing  make  the  setting  up  of  such  an  enormous  matrix 
( 332 2  elements)  quite  feasible,  although  certainly  not  trivial  when  one 
adds  to  it  the  effort  needed  to  obtain  the  values  of  [nm  from  which 
the  are  then  computed  by  Fourier  analysis  to  give  the  agm  according 

to  (2.13).  This  whole  operation  required  50  c.p.u.  minutes  in  the  AMDHAL 
470v/VI - 1 1  computer  at  O.S.U.,  using  double  precision  and  FORTRAN  H 
extended  in  the  highest  optimizing  mode.  Another  advantage  of  (2.51) 
is  that  it  permits  the  setting  up  of  the  normal  matrix  without  first 
having  to  form  the  observation  equations  matrix  A  ,  which  is  truly  gigantic 
(3.9  x  106  x  11  x  10**  elements). 

The  inversion  of  a  general  matrix  of  the  size  of  G  ,  since  the  number 
of  operations  required  to  invert  a  d  x  d  matrix  is  of  the  order  of 
d3  ,  would  be  impossibly  laborious,  requiring  something  like  (N+l)6  operations 
thereabouts.  As  only  the  small  non-zero  diagonal  blocks  have  to  be  inverted, 
and  only  half  of  those  are  different,  the  actual  number  of  operations  is 
of  the  order  of  ^(N+l)4  ,  or  some  96  (N+l)2  =  1.06  x  10  7  times  less 
than  for  a  general  matrix.  Using  the  same  computer,  compiler,  etc.  mentioned 


,  or 
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above,  the  inversion  of  G  required  some  15  min  of  central  processor  unit 
time.  This  is  still  a  very  large  number  of  calculations.  However,  because 
the  inversion  involves  the  processing  of  independent  blocks,  each  relatively 
small,  the  rounding  errors  due  to  finite  register  length  (64  bits  each  double 
word)  are  confined  so  they  can  not  affect  the  results  in  any  appreciable 
way.  Another  property  of  a  block-diagonal  matrix  is  that,  as  the  blocks 
are  created  and  inverted  independently,  the  whole  procedure  is  ideally 
suited  for  parallel  processing. 

The  programs  used  to  form  and  invert  G  are  documented  and  listed 
in  appendix  B  . 


2.10  The  Existence  of  G_1 


Differential  measurement-,  such  as  the  relative  line  of  sight  velocity 
between  the  satellites,  tend  to  be  associated  with  observation  and  normal 
matrices  that  are  rank -defficient,  so  it  is  reasonable  to  wonder  whether 
the  inverse  of  G  ,  so  important  to  an  error  analysis,  does  in  fact  exist 
at  all.  As  shownhere,  th^s  is  not  an  entirely  idle  question,  because 
there  are  cases  where  G  is  singular,  though  fortunately  not  with  the 
mission  parameters  chosen  in  this  study. 

Two  trivial  examples  of  singular  configurations  are  \p  =  0  ,  when 
thewhole  matrix  A  vanishes,  and  G  with  it,  and  ip  =  it  ,  when  all 
the  columns  of  A  corresponding  to  odd  harmonic  degrees  are  zero.  If 
A  is  singular,  G  is  singular  too,  so  G-1  does  not  exist  if 

Az  =  0 


for  some 

z  t  0 

Calling  ^nm  to  the  element  of  z  corresponding  to  in  c  ,  then, 

according  to  (2.13),  (2.22,a-b),  and  (2.24),  the  elements  y,  of  the 
vector  y  =  A^z  have  the  form 


N 

Y 


8=0  m=0  p=z 


I 


pm 


{ (pw  ±  mfi)ti 


(2. 52, a) 


where 


1  N 

i  i  «p 

a=Q  j =m  H 


f(Aa)  Z 


jm 


(2.52,b) 


with  1  sin(pu;  ±  mft)Aa 

f<4«>  *  (pST "  !i.cos(pu.*mSi)ia'  1 

ajm  =  h™  [(j+l)  cosp-|  sin-|-  +  p  sin  P  \  cos  (expression  (2.13)) 


Expression  (2. 52. a)  is  a  Fourier  expansion  with  coefficients  ,  and  it 
can  be  zero  only  if  all  such  coefficients  vanish,  as  long  as  the  sampling 
rate  is  higher  than  the  highest  frequency  in  the  expansion,  as  it  is 
the  case  here  according  to  paragraph  (2.7).  In  general,  some  elements 
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of  z  can  be  zero,  so  assume  that  Z'j*m  is  the  element  of  highest  degree 
and  order  that  is  not  zero.  In  such  case,  (nu,  +  ntf)  is  the  highest  angular 
frequency  in  y-j  .  The  only  terms  of  this  frequency  are  fp^  sin(nx  +  mft) 
and  f°nm  cos(nu)  +  mf>)  which  cannot  cancel  each  other  out.  For  both 
to  be  zero  for  all  ti  must  be  fpm.  *  ^nm  =  0  •  The  other  relevant 
factors  involved  in  (2.52,b)  are:  Z^m  ,  which  is  assumed  to  be  not  zero; 
n™  which  is  the  Fourier  coefficient  of  of  highest  frequency 

and  cannot,  for  this  reason,  be  zero;  sin(pu.  +  mf 2)  Aa  and  1  -  cos(n^  +  mH)Aa 
with  Aa  >  0  ,  which  cannot  be  simultaneously  zero  because  for  this 

to  occur  must  be  (n  m.C)At  i  (nu.  +  mfjAa  =  k2r ,  as  Aa  <  At  ,  but  this 
is  precluded  by  the  sampling  ratio  At~  being  more  than  twice  the 
highest  frequency  (nu+m.Q)/2Tr  and,  therefore  ,  (nu+m$Aa  <  (n^  +  rrp)  :t<  2m  . 
The  only  way  in  which  the  terms  of  frequency  (nu.  +  rrt")  can  vanish  for  all 
t-  ,  for  a  given  y  ,  is  that 

(n+1)  cos  n^sin^  +  nsinny  cos  ^  =  0  (2.53) 

for  some  n  in  2  <  n  N  ,  as  n  =  0  and  n  =  1  have  been  excluded 
from  the  model  (paragraphs  2.4  and  2.5).  This  is  a  necessary  condition 
for  G_1  not  to  exist.  Conversely, 

(n+I)  cosn^.  sin  i +n  sinn  ^  cos /  0  (2.54) 

for  all  n  in  2  s  n  *  N  is  a  sufficient  condition  for  G-1  to  exist, 
or  G  not  to  be  singular.  Fortunately  (2.54)  is  met  by  all  n  in 
2  <:  n  <.  331  with  the  chosen  inter-satellite  distance  of  300  km  (ip  approx. 

0.046  rad).  According  to  (2.53),  the  critical  values  of  ^  at  which 
(2.53)  is  met  for  some  n  in  the  band  must  be  isolated  points  in 
0  .<  0  i  2  .  If  ip  coincided  with  some  0  ,  and  this  were  truly  a  point 

of  singularity  (remember  that  2.53)  gives  only  a  necessary  condition), 
then  one  could  choose  ip  =  $  +  6y,  with  cy  a rb i tra r i ly  sma  1 1 ,  and  G*1 
would  still  exist.  From  a  practical  point  of  view,  G  should  become 
increasingly  ill-conditioned  as  5y  -*•  3  and  singularity  is  approached, 
so  eventually  it  would  be  impossible  to  calculate  G_1  numerically  even 
when,  theoretically,  the  inverse  does  exist.  For  that  reason  the  stability 
of  the  numerical  inversion  of  G  should  be  checked  by  obtaining  the  dif¬ 
ferences 

5z  =  z  -  G-'Gz  (2.55) 

where  G^1  is  the  computed  inverse  and  z  an  arbitrary  Nc  vector  whose 
components  are  chosen  from  a  sequence  of  random  rumbers.  The  stability  can 
be  judged  from  the  number  of  significant  figures  that  z  and  G^G  z  have 
in  common. 

For  small  ip  ,  as  it  is  the  case  with  intersatell ite  separations, 
the  first  term  in  (2.53)  can  be  much  smaller  than  the  second  term,  so 
the  necessary  condition  for  singularity  can  be  written,  after  some  obvious 
simplification 

sin  n  -jjf>  *  0  (2.56) 

which  is  met  for  those  critical  values  $  where  n<p  =  0  ,  2ir  ,  4-rr  ,  .  .  . 
k2*  ,  .  .  •  Looking  at  the  derivation  of  (2.13),  on  which  (2.53)  is  based, 
one  can  see  that  the  term  ignored  corresponds  to  the  effect  of  the  radial 

-34- 


component  of  the  Inertial  acceleration  on  the  line  of  sight  velocity.  In 
a  flat  Earth  this  component  would  be  always  perpendicular  to  the  line  of 
sight  and  have  no  effect  at  all,  so  (2.56)  could  be  applicable  to  a  flat 
Earth  geometry.  Breakwell  (1979)  has  used  a  flat  Earth  approach,  his 
results  showing  peaks  in  the  noise  to  signal  ratio,  or  percentage  error, 
at  those  spatial  wavelengths  y  satisfying  the  condition  ±  =  k >  with 
k  integer.  The  spatial  wavelength  on  the  sphere  for  a  harmonic  of  degree 
n  is  ;  replacing  y  with  -$■  in  "Breakwell 's  condition"  =  k-y  , 

one  gets  (2.56).  A  singularity  in  the  operator  involved,  the  G  matrix 
for  instance,  would  result  in  the  relative  error  being  infinite  at  some 
wavelength,  indicating  complete  loss  of  information  at  t.ie  corresponding 
frequencies  due  to  the  differencial  nature  of  the  measurements.  While 
such  singularities  were  not  found  in  the  spherical  Earth  analysis  reported 
here  because  of  the  parameters  chosen,  there  were  nevertheless,  some 
very  gentle  and  broad  ripples  in  the  relative  error  as  a  function  of 
the  harmonic  degree  n  ,  with  maxima  at  those  degrees  where  rty  was 
closest  to  k2n  ,  as  the  reader  can  see  by  looking  carefully  at  the  results 
listed  in  appendix  C  . 


2.11  Least  Squares  Collocation 

In  general,  a  linear  estimator  of  £  from  2  (obs )  has  the  ^orm 

i  a  F  ii2(0bs)  =  +  n)  (2-57) 

where  F  is  the  Nc  x  Np  estimator  matrix.  In  the  case  of  least  squares 
adjustment,  as  the  reader  may  remember,  this  matrix  was  called  Fp  (expres¬ 
sion  (2.41)).  As  the  least  squares  estimator  of  paragraph  2.8  is  unbiased 
when  the  model  X12  =  is  perfect,  the  error  in  c  is  due  purely 
to  the  data  errors: 


e  =  c  -  c  =  -Fn  (see  (2.42)) 


The  purpose  of  least  squares  adjustment  is  to  minimize  the  covariances 
of  the  elements  of  e^  or,  equivalently,  the  trace  of  the  error  matrix 


E 

n 


■  E  (e  eT  }  «  E  (Fe  el FT} =  FE  {e  e  } 

1 — n  — n  1  1  — n  — n  1  1-hi  — nJ 


Ft  =  F  D  FT 


(2.58) 


Not  all  linear  estimators  are  unbiased.  In  general 
«(,  8  c  ■  F  vJ2  /■  0  (2.59) 

where  e^  is  the  Nc  -  vector  of  bias  errors  in  £  .  The  gravity  field 
(£  ,  VjH  is  "deterministic"  (albeit  unknown),  while  the  noise  £  is 
"stocTiastic" ,  and  the  errors  in  £  that  one  and  the  other  give  rise 
to  are  also  "deterministic"  and  "stochastic",  respectively.  What  this 
difference  boils  down  to  is  that  the  noise  can  be  interpreted  mathematically 
as  a  random  process  while  the  field  and  the  bias  cannot.  To  obtain  estima¬ 
tors  that  minimize  the  total  error  (bias  plus  the  effect  of  the  noise) 
a  method  known  as  least  squares  collocation  has  been  developed  (Moritz,  1967,  Krarup 
1969)  that  treats  both  parts  of  the  error  in  a  way  that  is  formally  very 
similar,  by  using  the  operator  E(  }  for  the  noise,  and  the  operator 
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M  {  }  (average  over  all  rotations)  for  the  bias.  E  {  }  is  a  stochastic 
operator,  while  M  {  }  is  a  geometrical  operator.  If  is  the  bias 
in  the  estimate  and  er  the  effect  of  the  noise,  then  the  error 

measure  or  "covariance"  to  be  minimized  is 

*  M  K  }  +  E  } 

This  is  a  hybrid  error  measure  (geometric  plus  stochastic),  but  has  the 
advantage  of  being  quadratic  while  and  er  are  1  inear  functions 
of  vj2  and  £  ,  which  simplifies  the  mathematical  treatment.  The  rotations 
involved  in  the  application  of  M  {  }  are  those  of  the  system  of  coordinates 
and,  "attached"  to  it,  as  it  were,  of  the  pattern  in  which  the  data  are 
sampled  (i.e.,  the  orbits).  After  each  rotation  both  measurement  sites 
and  coordinates  are  different,  and  so  are  the  values  measured  and  the 
estimates,  together  with  the  bias.  M  {  }  averages  the  error  over  all 
such  "possible  outcomes"  (one  outcome  =  one  rotation).  If  the  statistics 
of  the  noise  £  are  gaussian,  so^are  those  of  cn  .  If  £5  =  0  ,  one 
can  determine  the  likelihood  of  C$m  being  within  so  many  "sigmas"  of 
the  true  value,  if  one  knows  only  the  covariance  and  the  mean  value  of 
the  data  errors.  With  the  "geometrical"  statistics  based  on  M  {  }  one 
would  need,  in  general,  not  only  M  {e^}  ,  but  the  higher  moments  M  {e&}  , 

M { e b }  ,  ...  as  well,  as  there  is  no  great  reason  at  present  to  believe 

that  the  gravity  field  is  sufficiently  "gaussian"  to  ignore  them.  Such 
moments  can  be  "propagated"  (in  the  sense  in  which  covariances  are  "propagated") 
from  the  corresponding  moments  of  the  geopotential.  The  second  moment 
is  the  covariance  function  (the  Legendre  transform  of  the  spectrum  a*) 
which  depends  only  on  the  geocentric  angle  between  two  points  on  the 
same  sphere.  Higher  moments  are  functions  of  many  such  angles.  Deter¬ 
mining  the  second  moment  or  covariance  from  empirical  data  is  not  easy, 
obtaining  the  higher  moments  must  be  even  harder,  so  probably  working 
with  these  moments  to  obtain  intervals  of  confidence  is  not  practicable, 
at  least  at  present.  In  practice,  the  choice  of  any  approximation  technique, 
such  as  collocation,  should  depend  on  how  well  it  works  for  the  sort 
of  problem  to  be  solved  with  it.  In  the  case  of  geopotential  coefficients 
determination,  this  author  has  conducted  several  numerical  tests  (Colombo, 

1981)  and  found  out  that  the  square  of  the  actual  errors  in  £  can  be 
very  close  to  their  hybrid  a  posteriori  covariances  when  estimating  £ 
with  collocation.  The  data  considered,  however,  were  point  and  mean 
gravity  anomalies,  not  SST  data.  For  a  discussion  of  the  theory  and 
applications  of  collocation  in  geodesy,  the  reader  could  see  Moritz  (1980). 

To  treat  biases  and  Dropagated  errors  in  a  way  that  is  formally  the 
same,  one  could  define  a  "variance-covariance"  matrix  for  the  biases 


As  already  mentioned  in  paragraph  1.2,  expressions  (1.15,a-c), 

"  tO  ■  0  ;  H  IICJ/l  ;  M  (C“m  CjM  -  0  if  o.  *  S  ,  n  f  k  ,  or 


m  f  q  ,  so  the  matrix 
C  =  M  {£  £T } 


(2.61) 
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is  diagonal.  Each  diagonal  element  in  C  equals 

Sap,  (2n+l)_1  if  n  i  M 
a2  (2n+l)_1  if  n  >  M 

where  the  Sa^  correspond  to  the  discrepancies  between  the  actual 

coefficients  and  those  of  the  reference  model  (paragraph  1.2).  Replacing 
(2.59)  in  (2.60)  and  multiplying  out  one  gets 

Eb  =  M  {(c  -  F  y12)  (c-Fvj^rMjcJ^c  v[2  FT)  +  M  {Fy12yL  FT} 

=  C  -  2CATFT  +  FACATFT  (2.62) 

because  vX2  =  Ac  .  The  error  measure  to  be  minimized  is  the  trace  of  the 
combined  error  matrix 

Et  =  F.  +  En  =  C  -  2CATF  +  FACAV,  FDFT  =  C-  2CATFT +  F (ACAT + D)  FT 
”  b  n  (2.63) 

according  to  (2.58)  and  (2.62),  or  tr  {Ej }  .  The  optimal  estimator  matrix 
minimizes  tr  {Ej}  ,  so  it  must  satisfy  the  matrix  equation 

*  ly-  {Ej}  =  F  (ACAT  +  D)  -  CAT  «  0  (null  matrix)  (2.64) 

This  equation  is  known  as  the  "normal"  equation.  The  solution  is 

F  =  CAT  (ACAT  +  D)'1  (2.65) 

Using  the  matrix  identity 

CAT  (ACAT  +  D)-1  -  (ATD_1A  +  C-1)  AT  D"1  (2.66) 

(see,  for  instance,  Uotila  (1967),  equation  (29))  the  optimal  estimator 
is,  finally,  according  to  (2.46),  (2.57),  (2.65),  and  (2.66), 

£  *  (atd-m  +  c-1)-1  AT  d-1  v12(obs) 

=  (G  +  C-1)  AV1  v12(obs)  (2.67) 

which,  except  for  the  term  C_1  ,  is  the  same  as  expression  (2.39) 

(with  P  *  D”  ;  for  leastsquares  adjustment.  If  all  the  degree  variances 
ojj  of  the  potential  are  non  zero  for  2  s  n  s  N  ,  C*1  exists,  G  +  C_1 
is  positive  definite  (sufficient  condition  for  F  in  (2.65)  to  minimize 
tr  {Ej}  ),  and  (G  +  C"1)"1  exists  as  well,  regardless  of  whether  G 
is  positive  definite  or  singular.  This  is  the  case  in  the  present  study, 
so  the  problem  discussed  in  paragraph  2.10  is  not  relevant,  at  least 
in  theory,  to  least  squares  collocation.  Ill -condi t ion i no  in  G  +  C_1» 
as  distinct  from  singularity,  is  another  matter.  Fortunately,  the  stability 
of  the  inversion  of  G  +  C"^  was  as  good  as  for  G  alone.  As  C'1  is 
diagonal,  G  +  C"1  retains  the  block-diagonal  structure  of  G  . 

Expression  (2.63)  can  be  transformed  as  follows,  when  F  satisfies 
(2.64): 


CaS  _ 

nmkq 
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(2.68) 


£T  =  C  -  FAC  =  C  -  (ATD"*A  +  C'1  pAVMC 
=  [I  -  (ATD_1  A  +  c-1)  aVa]C 
=  (aVm  +  C"1  )[aV1ATC"1  -atd_1a]c 
=  ( atd-1  A  +  C_1  )-L  I 
=  (G  +  C-1  P 


so  (G  +  C-1)"1  is  the  "hybrid"  variance-covariance  matrix  of  the  estimates, 
corresponding  to  G-1  in  least  squares ,  from  which  it  differs  only  in 
the  term  C-1  .  Expressions  such  as  (2.67)  and  (2.68)  have  been  used 
already  in  satellite  geodesy  for  modelling  the  geopotential  (Lerch  et  al., 

1977). 

In  general,  the  bias  matrix 

Eb  =  (ATD_1A  +  C-1)’1  -  =  (G  +  C'1)'1-  (G  +  C-1pG(G  +  r1  )_1 

=  (G  +  C-1  PC"1  (G  +  C'1  )_1  (2.69) 

is  not  zero,  so  the  optimal  estimator  of  collocation  can  be  a  biased  es¬ 
timator.  Before  the  introduction  of  collocation  in  geodesy,  the  deliberate 
use  of  biased  estimators  in  this  discipline  was  unusual.  In  other  fields, 
statistics  for  instance, biased  estimators  have  had  considerable  application,  for 
example  in  so-called  "ridge  regression"  (see  Bibby  and  Toutenberg,  1977) 
which  is  formally  very  similar  to  least  squares  collocation.  The  intentional 
use  of  biased  estimators  is  quite  common  in  communications ,  control  engineering 
and  in  data  processing  generally,  where  the  Wiener  and  Kalman  filters 
have  many  uses.  Both  types  of  filter  are,  in  fact,  biased  estimators 
whose  mathematical  formalism  resembles  that  of  collocation.  But  biased 
estimators  are  extensively  used  in  everyday  life  too:  every  radio  and 
audio  system  contains  filters  that  "estimate"  the  desired  component  in 
the  input  signal  (voice,  music)  by  rejecting  those  frequencies  at  which 
unwanted  signals,  or  noise,  are  most  prevalent,  and  reinforcing  those 
where  the  desired  ones  dominate  .  While  the  output  is  made  inteligible, 
the  signal  itself  has  been  distorted  because  those  frequencies  where 
it  overlaps  the  noise  have  been  smoothed  out  while  the  others  have  been 
boosted.  If  no  noise  were  present,  the  true  signal  would  not  be  identical 
to  the  "estimate"  or  audible  output,  the  difference  being  a  bias.  The 
example  is  quite  germane,  as  the  amplifiers  and  filters  in  question  are 
usually  linear,  like  the  estimators  of  collocation. 

Estimating  the  covariance  function  or  the  power  spectrum  of  the 
gravitational  potential  from  existing  data,  always  incomplete  and  imperfect, 
is  a  problem  that  entails  some  deep  and  difficult  theoretical  questions 
(Moritz,  1980).  From  a  practical  point  of  view,  one  could  estimate  the 
degree  variances  needed  for  C  in  (2.67)  and  (2.68)  from  the  SST  data 
itself  by  using  expressions  (1.18)  and  (1.19,a-c)  to  set  up  observation 
equations  relating  the  average  spectral  power  S™  (observable,  in  principle) 
to  the  unknown  oft  (0<ncN),  and  then  adjusting  the  latter  in  the 
manner  proposed  by  Wagner  and  Colombo  (1979)  for  high-low  SST  data.  The 
short-arc  spectral  technique  discussed  there  for  circular  orbits  was 
later  extended  by  Wagner  (1980)  to  elliptical  orbits  as  well. 
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In  the  last  analysis,  how  well  the  cr2n  must  be  known  depends  on  how 
sensitive  the  estimates  and  the  posteriori  "hybrid  covariances"  are  to 
errors  in  the  adopted  values  of  the  degree  variances.  In  section  3  there 
is  a  comparison  of  results  obtained  using  two  different  models  for  the 
power  spectrum,  showing  low  sensitivity  even  when  the  discrepancy  between 
the  spectra  is  considerable,  so  this  may  not  be  a  serious  problem. 


Throughout  the  discussion  one  has  considered  the  estimation  of  the 
actual  potential  coefficients  C^m  instead  of  the  scaled  .  The 
only  real  change  needed  is  in  expression  (2.61),  where  the  scaled  elements 
of  the  diagonal  of  C  are 


c 


aa 

nmnm 


(GM):  a2n  5an  ^2n+1)  1 
^0+*  x  a2n  (2nt-l)J 


n  <  M 
n  >  M 


(2.70) 


After  inverting  the  thus  modified  (G  +  C_i )  matrix,  the  actual  a  posteriori 
degree  variances 


*4*  l 


a=0  m=0 

can  be  obtained  from  the  scaled  a  posteriori 


,  where 
variances 


cta  =  e  -  ta  , 
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with  the  relationship 
,  R2n+4 

°n  =  iGMl^a211  °  n 


(2.71) 


2.12  Accuracy  of  the  Computed  Geoidal  Heights 

One  of  the  main  geophysical  applications  of  GRAVSAT  should  be  to 
provide  an  accurate  description  of  the  geoid  at  sea,  as  reference  surface 
for  oceanographic  studies  of  currents,  eddies,  etc.  The  geoid  height, 
or  difference  in  ellipsoidal  height  between  the  geoid  and  the  reference 
ellipsoid  ,  is 
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(2.72) 


assuming  that  all  masses  are  contained  inside  the  ellipsoid.  When  this 
is  not  the  case,  and  the  correspond  to  the  expansion  of  the  geopotential 

outside  the  bounding  spheren™i.e.  those  that  can  be  sensed  by  SST  data 
for  2  <  n  s  N),  expression  (2.72)  can  be  said  to  correspond  to  the  free 
air  geoid.  As  the  geoid  experiences  periodical  and  secular  variations, 
what  is  discussed  here  should  be  regarded  as  an  average  geoid  for  a  given 
period  of  time  (say,  during  the  GRAVSAT  mission),  but  this  average  can 
be  corrected  for  the  main  fluctuations,  such  as  tidal  effects,  to  obtain 
the  instantaneous  geoid. 


Using  the  estimated  potential  coefficients  one  can  compute  the  first 
(N  +  l)2-  3  terms  in  (2.72)  (the  band-limited  assumption  does  not  apply 
at  sea  level) 


»  (r  .♦,»)•  r  I  !  £  <♦.*>  £>" 


a=0  »-  2  m*0  m  ™ 
and  the  difference  between  (2.72)  and  (2.73)  is  the  error  in  geoid  height 


(2.73) 
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The  square  of  this  error  is 

«N(r,».*)‘.  r‘  l  l  l  l  l  l  CXf nmfkq(?2n- 

a=0  3=0  n=2  k.=2  m=Q  q=Q  nm  Kq  nm  *q 


where  = 
nm 
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If  the  Cnm  are  obtained  by  least  squares  adjustment  »  the  variance  of  the 
errors  would  be 


E  i<5 N ( r ,  <?,  A)2  j=o?<$N(r,  4> ,  A ) 

Keeping  r  constant  and  averaging  a2of<(r,4)  ,  A  )  over  the  unit  sphere 
(i.e.,  with  respect  to  <p  and  A  alone) 
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because  of  the  orthogonality  properties  of  spherical  harmonics.  In  an  entirely 
similar  way  one  can  arrive  to  exactly  the  same  expression  for  the  mean 
square  error  when  the  Cnm  are  estimated  by  collocation,  except 
that  at  hi  is  now  a  hybrid  degree  variance  in  the  sense  explained  in 
the  previous  paragraph.  As  the  error  depends  on  r  ,  and  r  is  not 
constant  on  the  ellipsoid,  one  could  settle  for  an  "average"  error  by 
choosing  r  *  a  ,  where  a  is  the  mean  Earth  radius,  so 
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*  a 
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m*2 


"ii 


+  l  o2 

n=N+l 


(2.76) 


which  is  the  expression  used  in  section  3.  Notice  that  the  correlations 
between  the  errors  cC^m  »  which  exist  when  they  have  the  same  order 
m  because  (G  +  C-1)*1  has  the  same  block-diagonal  structure  as  (G  +  C_1), 
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are  eliminated  from  (2.75)  and  (2.76)  by  the  orthogonality  of  the  f^m(<i>,lO  . 

To  obtain  the  accuracy  of  the  geoid  height  computed  with  the  SST  derived 
model,  one  has  to  know  the  a  posteriori  accuracies  of  the  coefficients, 
represented  by  the  cr2e^nm  ,  and  also  the  actual  degree 
variances  of  the  field  for  n  >  N  .  The  first  we  know,  to  some  extent, 
as  formal  accuracies  derived  from  the  diagonal  elements  of  (G  +  C"1)'1  ; 
the  second  are  largely  unknown,  as  they  correspond  to  the  part  of  the 
potential  not  analyzed  yet.  Therefore,  one  must  adopt  some  approximate 
model  for  this  spectral  "tail"  between  N  and  00  ,  baseJ  largely  on  indirect 
evidence,  which  means  that  the  estimated  contribution  of  this  "tail"  to  3N  , 
or  truncation  error,  must  depend  on  rather  weak  assumptions.  A  further 
problem,  of  more  theoretical  nature,  is  that  (2.72)  is  valid  only  when 
the  potential  can  be  expanded  in  spherical  harmonics.  When  there  are 
masses  above  the  geoid,  as  it  is  the  case  in  reality,  the  expansion  is 
known  to  converge  only  outside  the  smallest  sphere  bounding  all  those 
masses,  or  external  bounding  sphere,  and  not  necessarily  at  the  Earth's 
surface  itself,  where  it  is  actually  needed.  This  question  is  partially 
answered  by  a  theorem  first  proposed  by  Walsh  (1927),  reported  later  by 
Keldych  and  lavrentieff  (1939),  and  independently  stated  by  Krarup  (1969), 
according  to  which  there  are  always  spherical  harmonic  expansions  like 
the  one  in  (2.72)  that  approximate  uniformly  the  potential,  to  within 
any  arbitrary  accuracy,  on  and  above  the  Earth’s  surface.  A  consequence 
of  this  theorem  is  that  the  C$m  for  n  i  N  ,  with  N  finite,  in  the 
convergent  expansion  outside  the  external  bounding  sphere  (which  can  be 
detected  by  SST)  must  differ  from  those  in  some  of  the  internal  approx¬ 
imating  expansions  also  by  an  arbitrarily  small  amount.  In  this  sense, 
the  coefficients  detected  from  SST  data  can  be  considered  to  be  valid 
also  at  the  Earth's  surface. 


The  approximating  expansions  in  Walsh's  theorem  converge  uniformly 
to  the  functions  they  approximate  above  the  topography.  Their  analytical 
continuations  inside  the  Earth,  where  they  also  take  definite  values, 
have  no  physical  meaning  and  are,  thus,  free  from  the  natural  constraints 
that  determine  the  character  of  the  functions  represented  (potential, 
geoid  height,  the  estimation  error  9N  according  to  (2.74))  above  the 
topography  or  bathymetry  (near  the  sea  surface,  for  example).  The  free 
air  geoid,  say,  coincides  in  the  limit  with  the  true  geoid  whore,  tn-is  sur¬ 
face  runs  above  the  solid  Earth,  so  it  is  rather  smooth  and  gentle  in 
such  a  region,  but  it  can  be  argued  that  it  may  be  much  rougher  in  places 
elsewhere,  under  the  terrain.  The  mean  square  value  of  3N  in  (2.76) 
corresponds  to  an  average  over  the  sphere  of  radius  r  =  a  .  This  can 
be  interpreted  as  a  spherical  approximation  to  the  reference  ellipsoid 
after  all  masses  outside  it  have  been  removed  by  an  atmospheric  and  topo¬ 
graphic  correction.  This  is  the  interpretation  made  in  Section  3,  and 
it  involves  a  spherical  approximation  error  not  accounted  for  in  (2.76). 

One  can  also  regard  the  average  as  being  done  on  an  actual  sphere  that 
dips  in  and  out  of  the  topography  and,  more  importantly,  of  the  equatorial 
bulge,  so  a  correction  for  the  masses  above  this  sphere  is  too  large  to 
be  reliably  computed.  In  this  case  the  average  includes  the  squares  of 
true  errors  3N  above  the  Earth's  surface,  and  of  the  analytical  contin¬ 
uation  of  9N  below.  As  the  continuation  can  have  a  very  different  character 
from  that  of  3N  in  free  space,  one  may  wonder  just  how  meaningful  is 
this  average,  how  well  does  it  represent  the  real  errors,  particularly 
where  one  is  most  concerned:  near  the  sea  surface,  for  instance.  Finally, 
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the  presence  of  sharp  features  such  as  ridges,  trenches,  cordilleras,  etc., 
results  usually  in  strong  variations  of  the  geoid  which  may  cause  oscil¬ 
lations  in  the  truncated  series  of  the  estimated  N(r,  <f> ,  A)  of  the  type 
called  "Gibson's  phenomenon"  in  Fourier  series,  oscillations  that  may 
affect  the  practical  use  of  the  global  model  in  areas  that  are  among  the 
most  interesting  from  a  geophysical  point  of  view.  While  all  these  problems 
may  have  simple  solutions,  none  are  known  to  this  author  yet,  so  there 
may  be  some  need  for  improvement  in  our  present  understanding  of  these 
rather  basic  questions,  and  for  a  grain  of  salt  when  taking  some  of  the 
results  distilled  from  existing  theories.  With  such  reservations  in  mind, 
one  can  regard  the  undulation  estimate  of  (2.73)  as  optimal  in  the  sense 
that  if  the  SST  data  were  used  to  estimate  N  directly,  instead  of  finding 
the  Cftm  first,  the  result  would  be  exactly  the  same.  This  is  so  because 
tne  field  at  satellite  height  his  been  assumed  to  be  band  limited,  and 
in  such  case  any  linear  combination  of  the  Cftm  ,  such  as  N(r,  $ ,  X )  i  s 
according  to  (2.73),  can  be  optimally  estimated  by  replacing  the  Cftm 
with  their  optimal  estimates  in  the  linear  combination  (see  Colombo 

(1981),  paragraph  2.18).  The  truncated  part  (n  >  N)  is  regarded  here 
as  an  additional,  unknown  signal  on  which  there  is  no  information  in  the 
satellite  data,  while  the  terms  with  n  <  N  constitute  the  estimated 
variable. 


2.13  The  Effect  of  Some  Mission  Parameters  on  Coefficient  Accuracy 

The  diagonal  elements  of  the  inverse  of  the  normal  matrix,  o2e^nm 
in  the  sort  of  notation  used  for  the  elements  of  G  in  paragraph  2.9,  aepend 
on  the  values  of  certain  mission  parameters,  such  as  data  accuracy,  sampling 
interval,  mission  duration,  etc.,  according  to  simple  formulas.  Such 
formulas  enable  one  to  recalculate  the  accuracy  of  the  adjusted  coefficients 
for  different  values  of  those  parameters  without  having  to  repeat  the 
lengthy  computations  required  by  the  formation  and  inversion  of  the  normal 
matrix.  This  paragraph  contains  the  derivation  of  some  of  these  convenient 
formulas,  both  for  least  squares  adjustment  and  for  least  squares  collocation. 

(a)  Least  Squares  Adjustment: 

When  D  =  o2I  ,  where  I  is  the  Np  x  Np  unit  matrix,  the  normal 
matrix  is 

G  =  ATD_1A  =  o*2  ATI‘lA  =  o-2ATA 


and  the  variance-covariance  matrix  is 
G-i  =  cr2(ATA)”1 


so  the  diagonal  elements  of  G'1  are  related  to  o  by 
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(2. 77, a) 


The  elements  of  G  depend  on  the  sampling  interval  At  through  the 
number  of  samples  Np  =  r-£-  (expression  (2.51)).  Consequently,  the  elements 
of  G  are  inversely  proportional  to  At  ,  and  those  of  G*1  are  directly 
proportional  to  it,  so 
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(2.77 ,b) 


As  the  elements  of  G  are  directly  proportional  to  T  ,  the  length  of 
the  mission,  those  of  G”1  should  be  inversely  proportional  to  T  ,  provided 
that  T  =  k  T0  >  where  k  is  an  integer  and  T0  the  value  originally 
used  when  setting  up  G  (i.e.,  179  days  in  the  case  studied  in  Section 
3).  As  the  orbit  has  a  grand  period  T0  ,  extending  this  period  k  times 
results  in  the  orbit  being  repeated  also  k  times,  and  a  measurement 
at  any  given  location  being  taken  k  well,  so  the  adjustment 

contains  repeated  but  uncorrelated  measurements.  Consequently 


act 


=  (42-)a2e™  (T0) 

nmnm  '  T  nmnm 


(2.77,c) 


provided  T/f0  is  integer. 


For  small  changes  in  height  above  or  below  h  (the  height  used  to 
calculate  G  originally)  the  changes  in  satell  ite  angular  velocity  w  and, 
consequently,  in  the  shape  of  the  groundtrack  of  the  mid -point,  and  in 
the  distribution  of  the  measurements  along  it  are  small, too, and  can  be 
ignored.  If  \p  is  kept  constant  by  changing  the  separation  q  between 
satellites  so  sin-^-*  is  unchanged,  then  A  and  G  will 

vary  so  little  as  to  be  considered  independent  from  h  .  The  de-scaled 
diagonal  elements  of  G"1  are  then  related  to  h  as  follows 
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where  a  is  the  mean  Earth  radius.  Expressions  (2.77,a-d)  can  be  written 
together 
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in  accordance  with  (2.78),  provided  T/T0  is  integer,  |h-h0|  is  small 
compared  to  R  =  a+h  ,  and  \p  is  kept  constant. 


(b)  Least  Squares  Collocation: 

Rummel  et  al.  (1979)  have  shown  how  a  singular  value  decomposition 
of  the  normal  matrix  greatly  simplifies  the  recalculation  of  the  a  posteriori 
variances  with  different  levels  of  data  noise.  This  idea  can  be  adapted 
to  the  global  case  under  consideration.  It  follows  from  the  previous 
discussion  that 

G' (o,flt,T)  *  2u^4-  S’fro.AtB.To)  (2.80) 

where  G1 (oojAto.To)  is  the  normal  matrix  G  of  least  squares  adjustment 
pre-  and  post-multiplied  by  C*  and  calculated  with  the  parameter  values 
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Jo.  A t0 .  T0  .  If  A0  is  the  diagonal  matrix  whose  diagonal  elements  X. 


xi 


are  the  eigenvalues  of  G'(o0,  At0 ,  To),  if  all  the 
and  nonzero,  and  if  M0  is  the  matrix  whose  columns 
ponding  normalized  eigenvectors  of  G'(a0,  At0,  T0),  so'  m-j1  m-j  =  1 
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It  follows  that 
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The  eigenvalues  of  G‘(o,At,T)  +  I  are  X-j  +  1  ,  while  the  eigen¬ 

vectors  are  those  of  G'(a,At,T)  .  The  eigenvectors  of  a  real  symmetric 
matrix  with  distinct  eigenvalues  are  orthogonal  and  equal  in  number  to 
the  dimension  of  the  matrix,  so  M0  is  a  square  orthogonal  matrix  (the 
columns  are  orthogonal  unit  vectors,  if  one  assumes  that  all  the  X 1* 
are  different).  Consequently,  Mj1  =  MJ  ,  so 
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Therefore,  the  de-scaled  diagonal  elements  of  the  inverse  of  the  normal 
matrix, 
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(2.82) 

(2.83) 

(2.84) 


o'ehmnm  has  the  same  position  in  the  diagonal  of  (G  +  C*1)"1 
in  the  diagonal  of  A0  •  The  presence  of  C_1  in  the  expression 
of  the  normal  matrix  seems  to  preclude  a  simple  relationship  between 
the  a  posteriori  variance  and  h  like  that  of  (2.77,c),  because  C  depends 
on  R  =  a+h  ,  so  G'  =  C*  G  C*  and  A0  are  also  functions  of  h  (expressions 
(2.70)  and  (2.81). 

While  not  as  straightforward  as  (2.78),  (2.84)  is  relatively  easy 
to  apply,  compared  to  a  full  recalculation  of  the  normal  and  its  inverse. 

To  apply  this  expression  one  needs  to  have  the  eigenvalues  and  eigenvectors 
of  the  normal,  instead  of  its  inverse.  Obtaining  either  requires  much 
the  same  amount  of  computing. 


2.14  The  Right  Hand  Sides  of  the  Normals 

What  follows,  though  not  applicable  to  the  error  analysis  that  is 
the  main  subject  of  this  report,  is  relevant  to  the  question  of  data 
processing  discussed  in  Section  5. 

The  normal  equations  of  least  squares  adjustment  and  of  least  squares 
collocation  have  the  same  independent  terms  which,  in  vector  form,  c4n 
be  represented  by  the  Nc  -  vector 


(expressions  (2.47,b)  and  (2.67)) 
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An  individual  component  of  b  ,  corresponding  in  position  to  £. 
in  £  ,  can  be  written 

bnm  '  ^nm^  D  1 -12(obs)  =  ~  i|Q  anm^^  ^(obs) 

where  afjim  is  the  column  vector  in  A  corresponding  to  che  unknown 
(paragraph  2.7).  According  to  (2.24),  the  elements  of  a%m  are 
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The  C((poj  +  mf2)tj,Aa)  and  S((pw  -  mfi)t-j,Aa)  are  orthogonal  funcions 
of  ti  ,  according  to  (2.34,a-b),  so  the  data  vt^i)  .  can  be  expanded 
into  a  sum  of  these  functions: 

^12  (obs)  =  fQ  ^°pm  C((pl)  +  m«)t.^a)  +  v^  s((P-  +  mn)ti^a)  (2.87) 

P  m  r 

where  Nv  =  is_the  Nyquist  frequency  as  defined  in  paragraph  1.2,  while 
Np  =  i  uifl”  1Nj  Trhe  Vpm  >  Vpm  coefficients  can  be  obtained  from  the  or¬ 
dinary  Fourier  coefficients  in  Vi*H2,\  =  Ny  Nr  Aom  cos(pw  +  mft)ti  + 
i  Iods;  l  i  v" 
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''pm  =  Apm  (1-  cos(pu  +  mn)Aa)+  Bpm  sin(pu>  +  mS7)Aa 

According  to  (2.34,a-b),  (2.85),  and  (2.86,a-b),  the  elements  of  t)  are 
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(as  in  (2.20)) 


-45- 


where,  as  before,  the  upper  parts  of  the  curly  brackets  correspond  to  (n-m) 
even,  the  lower  to  (n-m)  odd,  w2  ±  =  f nm ±  mf A 2  ,  and  p  and  n  have 

the  same  parity.  ( 1-cos (p^ ±  mfi)Aa ) 

The  Fourier  coefficients  Cpm  ,  Spm  can  be  obtained  by  means  of  the 
Fast  Fourier  Transform  in  any  of  the  special  forms  designed  to  handle 
data  sets  too  large  to  store  in  the  central  memory  of  a  computer. 


2 . 15  Obi ique  Orbits 

The  error  analysis  reported  in  this  work  is  concerned  only  with 
polar  orbits,  which  provide  the  fullest  data  coverage  because  they  include 
the  polar  regions.  In  the  case  of  oblique  orbits,  the  inclination  of 
the  orbital  plane  with  respect  to  the  equator  is  neither  0°  nor  90°. 
With  oblique  orbits,  even  if  the  other  assumptions  in  paragraph  2.1 
are  maintained,  expression  (2.9)  for  the  line  of  sight  inertial  accel¬ 
eration,  which  is  the  starting  point  for  deriving  the  observation  equa¬ 
tions,  has  to  be  modified.  Calling  ,  as  in  Figure  2.2,  F  to  the  geo¬ 
centric  angle  between  the  ascending  node  and  a  point  P  along  the  orbit, 
and  L  to  the  longitude  of  that  node,  the  band-limited  gravitational 
field  potential  can  be  written  as  follows 
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where  the  upper  part  of  each  bracket  corresponds  to  (n-m)  even  and  the  lower 
part  to  (n-m)  odd;  i  is  the  inclination  angle  between  the  equator  and 

the  orbital  plane;  F  and  L  are  functions  of  time 
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e(t)  ~  -1.35  x  10'6  cos  i  rad  s'1  (2.91,c) 


is  the  rate_of  precession  of  the  node,  which  is  zero  for  polar  orbits 
( i  =  t-71  Fnmql’O  is  the  normalized  inclination  function 
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(2.92) 


where  k  is  the  integer  part  of  and  c  is  summed  over  all  values 

that  make  the  binomial  coefficients  nonzero.  For  the  derivation  of  (2.90) 
and  (2.92),  see,  for  instance,  Kaula  (1966),  Chapter  3  (the  notation 
is  somewhat  different). 


The  line  of  sight  between  the  satellites  lies  always  in  the  instan¬ 
taneous  orbital  plane.  The  inertial  acceleration  can  be  decomposed  into 
three  orthogonal  components:  the  radial,  the  normal  to  the  radial  in 
the  instantaneous  orbital  plane,  and  the  perpendicular  to  this  plane. 

Of  these,  only  the  first  two  accelerations  have  nonzero  projections  on 
the  direction  of  the  line  of  sight.  The  radial  acceleration  is 
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[{Snm}  cos((n-2q)F  + 
'‘nm 


*1) 


+ 


sin((n-2q)  F  +  mL)] 


(2. 93, a) 


and  the  normal  to  the  radial  acceleration  (in  the  instantaneous  orbital 
plane)  is 


i  3V  = 

r  Tf 


GM 

7~? 


N 

Jo 


l  <£)"  I  (n-2q)  C  sin((n-2q)F  +  mL) 

m=0  r  q=0  nmq  bnm 

+  {r™}  cos{ (n-2q)F  +  mL)]  (2.93,b) 

unm 


The  coordinate  L  is  common  to  both  satellites,  as  their  orbits  are 
coplanar.  Adopting  F  =  ,  the  midpoint  nodal  distance,  as  the 

other  independent  variable  ,r  =  R  is  fixed),  the  modified  inertial  line 
of  sight  acceleration  (inertial  accel .  minus  constant  term  due  to  even 
zonal s)  is 


a12(R,F,L)  «  -(|-7  (R.Mr.L)  +  -|J  (R,F+f,L))  sin|- 


1  (jt  (R,F — ^,L)  rc(R,F+ 


T 


L)  cos  -j  -  ac 


(2.94) 


Replacing  F  and  L  with  t  as  the  independent  variable  in  (2.93,a-b), 
substituting  the  Cftm  with  the  scaled  Cnm  »  introducing  p  =  n-2q  as 
subscript,  instead  of  q  ,  eliminating  all  terms  of  frequency  zero,  and 
replacing  the  result  in  (2.94),  one  can  finally  arrive  to  an  expression 
for  the  observation  equation  resembling  (2.24),  after  a  rather  laborious 
process  quite  similar  to  that  in  paragraph  2.2: 


Aa 


N 

Jo 
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l 

$ax(m,z) 


S((pco  +  mfi,)ti,Aa) 
a  ^-C((pw  + mfi' )ti  ,Aa) f 
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{S((po)- mQ' )t,,Aa); 
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C((pa)  +  rnfi')ti,Aa) 

J  1  ,  t  p1  +  m£>  ) 

‘S((Pu*nfi')ti,Aa)i 


(2.95) 


^  nni 
+  a  ,n± 


-C((pw- mfi1 ) t t  ,Aa )  ,  v 

(G|E)(i)  !  s((p^  -  m3'  1 1  i  ,Aa } r  (p“-mC'r  *  ♦  r* 


where 


(obs) 

v 


a°,(1Il|E)(t)  =  2pnm(d|E)(0  [(n+D  cosp-j-sin-^  +  p  sinp-jy-  cosy]  (2. 


96) 


While  their  coefficients  are  now  different  in  each  pair,  corresponding  pairs 
of  frequencies  p^  ±  n£  and  purnil*  ,  in  the  same  order,  appear  in  (2.24) 
as  they  do  in  (2.95).  As  a  result  of  this,  the  same  columns  in  A  are 
orthogonal  to  each  other  regardless  of  the  angle  ;  and,  therefore.the 
same  elements  in  G  =  c-: :A‘A  are  zero  as  before,  so  the  normal  matrices 
of  least  squares  adjustment  and  least  squares  collocation  have  the  block- 
diagonal  structure  shown  in  Figure  2.1,  paragraph  2.9,  for  oblique  as  well 
as  for  polar  orbits.  This  is  true  provided  that  w  and  TP  are  congruent, 
so  the  total  number  of  satellite  revolutions  Nr  and  the  total  number 
of  apparent  turns  of  the  Earth  in  a  node-fixed  system  of  coordinates, 
during  the  whole  mission,  Ng  ,  are  relative  primes. 

For  polar  orbits,  the  inclination  functions  have  the  property 


(p  t  0,  p  and  n  with  the  same  parity).  Comparing  (2.13)  and  (2.96) 
one  gets,  therefore, 

^  =  2Fnm(^)(-2-) 


(2.98) 


Consequently,  the  Fourier  coefficients  h™  of  the  Lnm  used  in  the 
case  of  polar  orbits  can  be  calculated  using  the  formulas 


/2n+T  (m=0) 

/2(2n-*-l )  (n-m) ! 
(n  +  m) ! 


min(k,^)  (2n  -  2t) !  . 

t~n  t!  (n-t)!  (n-m-2t)l  Z^n-t) 


*  l 


^n-m-2tj 


/  m 

((^rt_c) 


derived  from  (2.92)  with  \  -  y-  .  This  formula  is  an  alternative  to 
the  use  of  the  Fast  Fourier  Tr^hsform  as  explained  in  paragraph  2.7. 


(2.99) 


The  similitudes  between  (2.24)  and  (2.95)  indicate  that  the  computer 
programs  needed  to  carry  out  the  error  analysis  of  a  mission  where  the 
orbital  plane  is  inclined  with  respect  to  the  equator  are  very  similar 
to  those  for  polar  orbits  explained  and  listed  in  Appendix  B.  Reasoning 
once  more  as  in  paragraph  2.9.  one  arrives  to  an  expression  for  the  general 
element  of  the  G  matrix 
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3.  Numerical  Results 


This  section  presents  the  results  of  the  error  analysis  whose  theory 
has  been  given  in  sections  1  and  2  .  The  calculations  have  been  done 
in  accordance  to  least  squares  adjustment  and  to  least  squares  collocation, 
and  both  sets  of  results  are  shown  for  comparison. 


3.1  Spectral  Model  and  Error  Formulas 


The  degree  variance  of  the  error  is,  according  to  paragraph  2.13, 
n  1 

2  -  r  r 

°  c-n  L  L  ^nmnm 

m=u  j/=0 

~  -  "0. 

where  0  >,mnm  is  the  variance  of  the  estimated  coefficient  Cnm  •  The 

relative  error  per  degree  is 


(3.1) 


P 

n 


j2fn 

Jn 


where 


a* 

n 


n 


I 

m=Q 


(3.2) 

(3.3) 


is  the  degree  variance  of  the  n  harmonic  of  the  potential.  The  set 
of  all  for  0  <  n  <  °°  is  the  power  spectrum  of  the  potential.  This 
relative  error,  multiplied  by  100,  is  listed  as  percentage  error  per  degree 
in  the  various  tables  shown  in  this  section.  To  calculate  (3.2)  it  is 
necessary  to  know  the  power  spectrum,  the  a2  •  The  spectrum,  like  the 
potential  itself ,  is  not  entirely  known,  but  there  are  estimates  of  the 
o2  obtained  from  the  analysis  of  terrestrial  and  satel 1  ite  data .  For 
low  degree  harmonics  the  ^  can  be  calculated  using  (3.3)  with  the 
values  for  the  taken  from  one  of  the  existing  spherical  harmonic 

models  of  the  potential.  For  higher  terms,  the  "tail"  beyond  the  highest 
degree  whose  C^m  are  known,  one  can  choose  among  the  several  formulas 
for  oi  as  a  function  of  n  that  are  available,  each  based  on  a  different 
type  of  approximation,  or  on  different  data.  For  this  study  the  degree 
variances  up  to  degree  n  =  100  were  taken  from  a  spherical  harmonics 
model  complete  to  n  =  180  ,  obtained  by  Rapp  and  associates  at  OSU  from 
the  analysis  of  a  global  set  of  mean  l°x  1°  gravity  anomalies  using  quadrature 
formulas.  The  anomalies  themselves  were  obtained  from  gravimetry  and  altimetry 
by  means  of  least  squarescollocation.  As  a  further  step  the  anomalies 
were  adjusted  in  a  combination  solution  that  included  the  coefficients 
of  GEM-9  as  data.  For  a  report  on  this  adjusted  data  set,  see  Rapp  (1978) 
and  also  Rapp  (1979a).  GEM-9  is  described  in  (Lerch  et  al.,  1977).  If 
one  ignores  orbit  errors,  then  expression  (1.23)  is  reduced  to  an  identity 
between  the  time  derivative  of  the  residual  line  of  sight  velocity  and 
the  residual  line  of  sight  inertial  acceleration.  The  residual  inertial 
acceleration  corresponds  to  the  difference .hetween  the  true  field  and  the  field 
model  used  to  compute  the  orbit.  If  C$m  (Model)  a  coefficient  of 
the  model,  and  NM  the  maximum  degree  in  that  model,  then  the  geopotential 
coefficients  corresponding  to  the  residual  accelerations  are  =  Ctfm  - 

Crfm(Model )  f0r  n  <  NM  ,  and  C«m  for  NM  <  n  <  00  .  The  degree  valances 
of  the  power  spectrum  are,  therefore,  =  m ^  (cnm  *  Cnm™°°e1')2  for 

n  <  NM  ,  and  ^  (c*m) 2  for  NM  <  n  .  T?!ea  d2  are  the  degree 
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variances  corresponding  to  the  errors  in  the  reference  field.  As  these 
errors  are  not  known  exactly,  one  must  use  some  "likely  numbers"  instead, 
such  as  the  formal  error  variances  of  the  coefficients  obtained  during 
the  adjustment  of  the  reference  model.  This  criterion  has  been  adopted 
here,  the  error  variances  be.ng  those  of  the  first  NM  degrees  in  Rapp's 
model,  with  NM  =  20  .  For  n  >  NM  ,  the  as  implied  by  Rapp's  coef¬ 

ficients  have  been  used  up  to  degree  100  .  For  n  >  100  the  of,  have 
been  calculated  with  the  two-term  formula 

o2  =  g^p(n-  l)'1[Al(r.  +  A)'1Sln+2  +  A2(n  +  B)_1(n- 2)_1S2n+2]  (3.4) 

where  a  =  6371  km,  ^.  =  982026.41  mgal  A=1  ,  8  =  2  ,  Al  =  3.4050, 

A2  =  140.03,  SI  =  0.998006,  and  S2  =  0.914232  .  These  parameter 
values  were  obtained  by  Rapp  (1979b)  by  fitting  the  formula  to  the  empirical 
variances  of  his  180,  180  model  and  other  data.  In  summary: 

For  2  <n  a  NM  ,  the  error  degree  variances  f2nm  of  Rapp's  model; 

for  NM  <  n  <  100  ,  the  degree  variances  of  Rapp's  model; 

for  100  <  n  s  2000  ,  the  o ^  according  to  (3.4).  The  harmonic  con¬ 
tent  of  the  geopotential,  according  to  that  formula,  is  negligible  for 
n  >  2000  . 

The  error  in  geoid  undulation  due  to  the  errors  in  the  coefficients 
up  to  degree  n  has  been  calculated  with  the  formula 

a2eN  =  a2  T  c2et  (3.5) 

k=2  K 

If  no  coefficients  above  n  are  estimated,  the  total  error  must  be 

2000 

a2eNT  =cr2eN  +  a2  £  cr2.  (3.6) 

1  n  k=n+l  k 

(Compare  to  (2.76)). 

Expression  (3.6)  depends  on  the  tail"  of  the  spectrum  containing 
the  high  frequency  terms.  The  higher  n  ,  the  least  that  is  known  about 
on  ,  so  this  "tail"  is  the  least  reliable  part  of  the  spectrum  model, 
and  the  total  error  calculated  with  (3.6)  is  the  least  credible  among  the 
results.  Originally,  when  the  least  squares  collocation  results  were 
obtained,  the  total  error  was  not  computed,  though  its  calculation  was 
added  in  the  main  program  afterwards,  when  other  results  were  found. 

For  completness 1  sake,  the  punched  output  of  that  first  run  was  read  by 
an  auxiliary  program,  which  then  produced  the  printout  shown  in  Table 
3.2,  including  the  total  error  in  the  last  column,  and  also  the  full  listing 
of  Table  C.2  in  Appendix  C.A  minor  mistake  in  the  auxiliary  program  resulted 
in  values  of  the  total  error  that  are  incorrect.  To  obtain  the  "true" 
values  en  ,  the  listed  values  e'n  should  be  corrected  according  to 
the  formulas 

en  =  (e'2  +  0.1074)*  (3.7) 
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3.2  Results  According  to  Least  Squares  Adjustment 


Table  3.1  shows  the  accuracies  of  potential  coefficients  and  geoidal 
undulations  estimated  from  SST  data  collected  during  a  mission  whose  para¬ 
meters  were: 

Circular,  polar  orbit, 

satellite  height:  160  k.  «..i, 

intersatellite  separation:  300  km, 

accuracy  of  the  data:  /?  x  I0~5m  s'1 , 

averaging  interval:  4  s  , 

sampling  interval:  4  s  , 

length  of  the  mission:  179  days  , 

maximum  degree  and  order  in  reference  model:  20  . 

The  error  analysis  was  carried  out  according  to  least  squares  adjustment 
theory,  as  put  forward  in  paragraph  2.8.  The  first  column  shows  the  relative 
percentage  error,  which  is  the  ratio  defined  by  expression  (3.2)  multiplied 
by  100  .  Notice  that  above  n  =  270  the  error  exceeds  100%  consistently. 

The  size  of  the  the  errors  oecome  so  large  that  the  total  undulation 
error,  which  decreases  steadily  up  to  n  =  270,  according  to  the  last 
column,  begins  to  increase  quite  perceptibly  once  more.  The  last  but  one 
column  shows  the  error  up  to  n  ,  according  to  expression  (3.5).  The 
second  column  contains  the  values  of  -§$-2(-|-)2nc2en  *  the  variance  per 
degree  of  the  error  in  the  scaled  coefficients  (see  expression  (2.70) 
in  par.  2.11).  Notice  the  very  low  percentage  errors  for  n  <  100  ,  which 
are  much  the  same  as  those  shown  in  the  next  paragraph  for  least  squares 
collocation. 

As  explained  in  the  last  part  of  paragraph  2.10,  the  error  in  the 
recovered  coefficients  may  show  local  peaks  at  those  degrees  satisfying 
the  condition 

n^  =  2  m  k  ( k  =  1 ,  2 ,  .  .  .  ) 

Since  v  =  2sin_1(^)  ,  for  a  separation  of  300  km  the  maxima  should  occur 
within  the  band  0s  n  s  331  at  degrees  n  =  136  and  n  =  273  ,  approx¬ 
imately,  the  next  peak  above  the  band  being  at  n  =  410  .  The  listing 
in  table  3.1  is  too  coarsely  spaced  in  n  to  show  these  peaks,  which 
are  rather  narrow  features  .each  atop  a  broader  rise  that  surrounds  it, 
but  they  are  quite  clear  in  the  detailed  listing  of  Table  C.l  in  Appendix C. 

The  values  that  appear  both  in  Tables  3.1  and  in  Table  C.l  can  be 
modified  according  to  expression  (2.79)  in  paragraph  2.13,  to  obtain  the 
results  corresponding  to  missions  with  different  parameters.  Values  cor¬ 
responding  to  undulation  errors  (last  two  columns)  are  in  meters. 
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3.3  Results  According  to  Least  Squares  Collocation 

Table  3.2  shows  the  results  corresponding  to  least  squares  collocation, 
the  principle  of  which  has  been  explained  in  paragraph  2.11.  The  mission 
parameters  are  the  same  as  for  Table  3.1.  Comparing  both  tables  one  can 
see  that  the  collocation  accuracies  are  consistently  better  than  those 
for  least  squares  adjustment,  particularly  above  degree  n  =  210  .  The 
results  listed  in  the  last  column  of  this  table,  and  also  of  Table  C.2, 
Appendix  C  ,  have  to  be  corrected  according  to  expression  (3.7),  at  the 
end  of  paragraph  3.1.  When  contrasting  tables  3.1  and  3.2  one  must  bear 
in  mind  that  the  first  corresponds  purely  to  propagated  noise,  while  the 
second  contains  hybrid  variances:  noise  plus  bias.  The  peaks  in  the 
error  mentioned  in  the  previous  paragraph  are  much  less  noticeable,  in 
fact  the  one  at  n  =  273  has  disappeared  altogether  (Table  C.2),  and  only 
an  "acceleration"  in  the  increase  of  the  error  with  n  remains  in  the 

neighborhood  of  the  missing  peak.  Nowhere  the  error  exceeds  100%,  also 

according  to  the  theory.  1 

Table  3.3  lists  the  accuracies  corresponding  to  a  mission  with  the 
same  parameters  as  in  the  preceeding  tables,  except  that  the  accuracy 
of  the  data  is  now  four  times  worse,  i.e.,  /J  x  4  x  10~6ms_1  . 

Table  3.4  corresponds  to  the  same  parameters  as  in  Table  3.1,  but 

the  height  has  been  changed  to  220  km  above  the  Earth. 

Table  3.5  is  for  the  same  parameters  as  Table  3.4,  with  the  data 
noise  increased  to  v2  x  4  x  10-6m  s-1  .  The  results  for  Tables  3.3  through 

3.5  are  clearly  worse  than  for  Table  3.2,  as  can  be  expected,  because 
the  data  is  noisier,  the  signal  weaker  (higher  altitude),  or  both,  com¬ 
pared  to  the  case  of  Table  3.2. 

Table  3.6  corresponds  to  the  accuracies  that  would  be  achieved  if 
all  the  coefficients  in  a  given  degree  could  be  estimated  with  the  same 
accuracy  as  the  zonal  harmonic.  The  results  were  computed  with  a  modified 
version  of  the  main  program  of  Appendix  B  that  sets  up  and  inverts  only 
that  part  of  the  normal  matrix  corresponding  to  the  zonal s.  In  this  way, 
an  approximate  analysis  can  be  carried  out  at  much  less  cost  than  the 
complete  studies  of  tables  3.1  through  3.5.  The  mission  parameters  in 

3.6  are  the  same  as  in  Tables  3.1  and  3.2.  In  Table  3.7  the  separation 
between  satellites  has  been  changed  to  150  m  ,  and  in  Table  3.8,  to  600 
km.  The  results  in  Table  3.7  are  clearly  much  worse  than  for  Table  3.6, 
showing  a  particular  deterioration  at  low  degrees  due  to  the  differential 
nature  of  the  SST  measurements  that  tends  to  eliminate  the  low  frequency 
information,  so  the  percentage  of  the  error  in  the  estimated  coefficients 
increases  as  the  satellites  become  closer.  Table  3.8  shows  the  smallest 
band  error  and  the  best  accuracies,  but  these  are  quite  irregular  because 
the  error  peaks  in  this  case  are  spaced  closer  than  in  the  other  cases, 
and  are  much  more  prominent,  so  that  they  can  be  noticed  even  with  the 
coarse  spacing  in  n  used  in  the  table. 

From  the  results  given  in  this  paragraph  one  can  conclude  that  the 
orbits  should  be  as  low  as  possible,  and  the  separation  between  satellites 
as  wide  as  allowed  by  natural  limitations,  the  main  among  which  is  the 
need  to  keep  the  radar  beam  from  entering  too  deep  into  the  upper  layers 
of  the  atmosphere  at  the  mid  point.  Data  noise,  of  course,  should  be 
as  low  as  technically  feasible. 

^^The  optimal  estimator  should  not  be  worse  than  a  null  estimator 
(one  that  predicts  only  zeroes),  whose  error  is  always  100%. 
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Table  3.3  Parameters:  height  160  km,  separation  300  km,  o  =  /2  x 
Procedure:  least  squares  collocation 
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PERCENTAGE  ERROR  ERROR  VARIANCE< POT. >  HAND  ERROR  (  UNDUL. )  (  M)  TOTAL  UNDUL.  ERROR  <M) 
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Finally,  the  results  listed  in  Table  3.9  below  show  the  sensitivity 
of  the  relative  accuracies  of  the  estimated  coefficients  to  the  choice 
of  power  spectrum  model .  The  degree  variances  - ^  of  subroutine  NVAR, 
described  in  paragraph  3.1  and  in  appendix  B  ,  have  been  replaced  by  those 
obtained  according  to  a  two-term  model  of  the  form  described  by  expression 
(3.4),  but  with  the  following  parameters'  values: 

A  =  100  ,  B  =  20  ,  Hi  =  18.3906  ,  A2  =  658.6132  ,  51  =  0.9943667  , 

S2  =  0.S08949  . 

These  values  correspond  to  a  model  obtained  by  Jekel i  ("21" 
model.  Report  No.  275,  Dept,  of  Geodetic  Science,  The  Ohio  State  University, 
Columbus,  Ohio,  1978).  The  values  in  the  first  column  are  the  same  as 
in  Table  3.6,  those  in  the  second  column  correspond  to  the  "2L"  spectrum. 

As  in  the  case  of  Tables  3.6  to  3.8,  these  values  correspond  to  the  zona Is 
only. 


Table  3.9 


Relative 

Accuracies  with  Two  Different 

Spectral  Models 

n 

As  in  Table  3.6  {%) 

Jekel i's  "2L"  (! 

2 

.54171  x  10-1 

.55870  x  10“ 

20 

.74292  x  10‘2 

.74343  x  10- 

40 

.45779 

.36486 

60 

.10564  x  10-1 

.87091 

80 

.30675 

.20738  x  10- 

100 

. 10694 

.72841 

120 

.51907 

.382197 

140 

2.9581 

2.1529 

160 

1.1748 

.86360 

180 

1.4831 

1.1068 

200 

2.5173 

1.9143 

220 

5.2505 

4.0815 

240 

13.900 

11.106 

260 

53.941 

46.794 

280 

86.313 

82.380 

300 

70.181 

65.245 

320 

76.689 

73.190 

330 

76.327 

79.894 

The  dn  up  to  n  =  NM  =  20  are  the  same  for  both  columns  (subrou¬ 
tine  MODEL  in  appendix  B).  Clearly,  while  the  change  in  spectrum  does  make 
a  difference,  the  changes  have  little  influence  on  the  calculated  accuracies. 
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3.4  Accuracies  of  Different  Harmonics  of  the  Same  Deqree 


Having  come  across  the  argument  that  SST  data  collected  in  a  polar 
orbit  should  be  most  sensitive  to  the  zonals,  as  all  their  variation  occurs 
in  the  N-S  direction,  and  least  sensitive  to  the  sectorials,  with  the  tesserals 
falling  somewhere  in  between,  and  the  difference  from  zonals  to  sectorials 
being  quite  large,  the  author  included  in  the  main  program  statements 
to  list  the  accuracies  for  the  cosine  terms  (a  =  0)  of  all  potential  coef¬ 
ficients  of  degree  2  -s.  n  <  40  .  The  accuracies  for  n  =  30  are  listed 
below.  The  mission  parameters  are  those  for  Table  3.1,  and  the  principle 
used  is  that  of  least  squares  adjustment.  The  accuracies  correspond  to 
dimensionless,  scaled  coefficients  (par.  2.3). 


Table  3. 10 

Accuracies  of  potential  coefficients  of  degree  n  =  30 


m  ^nm 


0 

2.08  x 

10'13 

1 

1.46 

it 

2 

1.50 

It 

3 

1.51 

H 

4 

1.53 

n 

5 

1.58 

!• 

6 

1.61 

ll 

7 

1.68 

ii 

8 

1.71 

ll 

9 

1.78 

ll 

10 

1.83 

ll 

11 

1.90 

ll 

12 

1.97 

li 

13 

2.03 

ll 

14 

2.11 

ll 

15 

2.17 

ii 

16 

2.27 

il 

17 

2.33 

(l 

18 

2.44 

ll 

19 

2.49 

ll 

20 

2.62 

il 

21 

2.67 

ll 

22 

2.72 

ll 

23 

2.87 

:i 

24 

3.05 

n 

25 

3.09 

ii 

26 

3.31 

ii 

27 

3.35 

n 

28 

3.60 

ii 

29 

3.68 

ii 

30 

2.82 

#i 

The  values  above  are  typical  of  those  listed  for  other  degrees  in 
the  interval  2  s  n  $  40  .  While  there  are  fluctuations,  the  sensitivity 
of  the  adjustment  to  the  various  harmonic  terms  of  degree  n  =  30  does 
not  change  very  much. 


4.  Validity  of  the  Results 


The  numerical  results  of  the  previous  sectior  have  been  obtained 
under  the  simplifying  assumptions  of  paragraph  2.1.  Among  those  assumptions, 
numbers  (1)  to  (4)  define  an  orbital  geometry  more  regular  than  what 
can  be  found  in  reality,  and  assumption  (11)  disregards  all  sources  of 
error  other  than  SST  data  errors.  Of  the  remaining  ones,  (9)  and  (10) 
have  been  explained  already  in  the  first  section  of  this  report,  while 
(5)  to  (8)  merely  describe  what  an  ideally  successful  mission  would  produce 
in  terms  of  data.  Assumption  (3)  defines  an  idealized  Earth  rotation, 
where  there  are  no  fluctuations  in  the  angular  velocity  ft  due  to  such 
causes  as  tidal  friction,  redistribution  of  atmospheric  masses,  lack 
of  rigidity  of  the  Earth,  etc.,  and  no  changes  in  the  inertial  orientation 
of  the  spin  axis  due  to  precession-nutation  and  polar  wandering.  Changes 
in  .n  are  of  two  kinds:  short  term,  due  to  solid  Earth-atmosphere  inter¬ 
action,  etc  ..altogether  probably  too  small  to  matter,  and  secular,  due 
to  tidal  effects  caused  by  the  Sun  and  the  Moon,  also  very  small  and 
quite  predictable  from  long  records  of  observations.  Polar  motion  is 
also  well  modelled  from  very  long  series  of  observations.  The  main  effect 
of  changes  in  Earth  rotation  is  in  the  calculation  of  the  satellites' 
orbits,  as  orbital  errors  affect  also  the  accuracy  of  the  results.  Because 
of  the  existence  of  good  models,  this  effect  is  probably  negligible  in 
the  present  context.  So  this  section  is  going  to  consider  only  the  conse¬ 
quences  of  assumptions  (1),  (2),  (4),  and  (11)  on  the  credibility  of 
the  results  of  the  error  analysis.  The  basic  argument  is  that  the  actual 
data  set,  with  its  complex  three-dimensional  distribution,  can  be  reduced 
or  transformed  into  another  with  the  geometry  implied  by  the  assumptions 
and  with  the  same  information  content  as  the  original.  The  analysis 
of  this  transformed  data  set  can  be  done,  then,  in  the  manner  explained 
in  section  2,  the  accuracy  of  the  estimated  coefficients  being  much  the 
same  as  that  shown  in  the  previous  section. 


4.1  The  Geometry  of  the  Real  Orbit 

The  departure  of  the  gravitational  field  from  that  of  a  central 
point  masscauses  the  orbit  to  take  a  shape  that  is  not  exactly  circular, 
however  carefully  the  satellites  are  manouvered  into  it.  Most  of  the 
departure  from  cicularity  is  due  to  the  part  of  the  anomalous  field  that 
is  already  known  from  the  study  of  terrestrial  and  spacecraft  data,  and 
the  major  portion  of  this  known  departure  is  caused  by  the  second  zonal, 
which  is  almost  three  orders  of  magnitude  larger  than  any  other  harmonic. 
This  zonal  represents  most  of  the  effect  of  the  Earth's  equatorial  bulge 
on  the  geopotential.  The  result  is  a  wavy  motion  of  the  satellites, 
which  alternatively  run  above  and  below  the  meansphere  of  radius  R  (average 
orbital  radius),  gain  or  loose  ground  in  the  along-track  direction  with 
respect  to  a  perfectly  uniform  circular  motion,  and  also  move  with  respect 
to  each  other  because  of  their  different  positions  along  the  (more  or 
less)  common  orbit,  so  the  orientation  of  the  line  of  sight  is  not  always 
perpendicular  to  the  radial  direction,  but  fluctuates  about,  and  there  are 
also  variations  in  the  intersatell ite  distance.  The  non-zonal  terms 
of  the  harmonic  expansion  introduce  further  irregularities,  particularly 
in  the  across-track  direction,  so  the  orbit  does  not  lie  in  any  given 
plane  except  approximately.  The  effect  of  the  non-zonals  is,  however, 
of  the  order  of  meters,  while  that  of  the  second  zonal  amounts  to  several 
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kilometers.  The  discussion  that  follows  suggests  that  errors  of  a  few 
meters  have  negligible  consequence,  so  only  the  departures  due  to  oolatness 
need  to  be  considered.  If  all  measurements  had  been  taken  on  a  sphere  of 
radius  R  ,  the  line  of  sight  being  always  perpendicular  to  the  radial  direction, 
and  the  separation  between  satellites  constant,  but  the  actual  positions 
of  the  spacecrafts  had  been  displaced  in  latitude  and  longitude  from  the 
ideal,  periodical,  regularly  spaced  pattern  implied  by  the  assumptions,  one 
could  transform  this  "semi -perfect"  set  of  observations  into  a  "perfect" 
set  by  interpolating  horizontally  the  data  from  their  actual  positions  to 
their  ideal  ones.  As  far  as  the  signal  content  is  concerned,  this  inter¬ 
polation  can  be  exact,  because  the  information  is  band-limited,  in  the  sense 
of  paragraph  1.2,  and  the  sampling  is  very  dense.  That  such  interpolation, 
to  be  exact,  would  probably  require  the  use  of  all  data  values  to  create  each 
interpolated  value  is  of  no  consequence  here,  as  this  discussion  is  concerned 
only  with  the  the  existence  of  valid  transformation  procedures,  regardless  of 
whether  they  are  practical  or  not.  The  practical  aspects  of  the  matter 
are  left  for  section  5.  The  only  problem  would  be  that  the  interpolated 
noise  would  not  be  exactly  uncorrelated  and  of  constant  variance,  its 
departure  from  those  ideal  qualities  depending  on  how  much  the  true 
groundtrack  departs  from  the  ideal  one.  Here  one  can  only  assume  that, 
if  this  departure  is  not  too  large,  neither  would  be  the  change  in  the 
nature  of  the  noise  too  large.  At  this  stage  of  the  argument  it  will 
be  assumed  that  the  actual  positions  of  the  satellites  can  be  known  with 
negligible  errors  from  orbital  calculations,  the  effect  of  orbital  error 
being  treated  later  on  in  paragraph  4.3,  where  assumption  (11)  is  discussed. 

As  explained  there,  the  model  adopted  for  the  line  of  sight  velocity 
implies  that  such  errors  can  have  a  very  small  effect  on  the  estimated 
coefficients  as  long  as  they  do  not  exceed  a  few  meters,  an  accuracy 
achievable  with  existina  ^rbit  determination  techniques. 

The  transformation  of  the  actual  set  of  observations  into  a  set 
of  pseudo-observations  lying  on  the  mean  sphere  along  the  ideal  orbit 
can  be  carried  out  in  two  steps: 

(a)  a  vertical  reduction  of  the  original  observations  to  the  mean  sphere, 
where  an  intermediate  set  consisting  of  pseudo-obs.  with  the  following 
characteristics  is  created:  the  "midpoint"  between  "satellites"  lies  directly 
below  the  true  midpoint,  the  "line  of  sight"  is  oriented  North-South 

and  perpendicular  to  the  radial  direction,  the  "satellites"  are  at  the 
same  distance  from  each  other  in  all  the  pseudo  obs;  in  other  words: 
the  intermediate  pseudo  obs.  are  identical  to  ideal  observations  in  every¬ 
thing  except  their  arrangement  on  the  sphere; 

(b)  a  horizontal  interpolation  using  the  intermediate  data  set  to  form 
the  final  set  of  pseudo-obs.  at  regularly  spaced  points  along  the  ideal 
orbit. 

As  already  explained,  the  horizontal  interpolation  can  be  carried 
out,  in  principle,  exactly,  because  of  the  band-limited  nature  of  the 
gravitational  signal.  The  main  problem,  therefore,  is  the  vertical  reduction 
or  step  (a).  This  step  involves  a  downward  or  upward  continuation  of 


the  signal,  depending  <.n  where  the  satellites  happen  to  be  with  respect 
to  the  sphere,  and  a  correction  for  the  change  in  their  relative  position, 
which  results  in  a  variable  intersatell ite  distance  and  line  of  sight 
direction.  Both  the  continuation  and  this  correction  can  be  considered 
together,  as  an  overall  operation. 

4.2  Vertical  Reduction  to  the  Mean  Sphere 

The  discussion  can  be  simplified  by  considering  observations  of 
line  of  sight  relative  acceleration,  ai2  ,  rather  than  line  of  sight 
velocity,  v12  .  The  signal  being  band-limited,  it  is  possible  to  dif¬ 
ferentiate  the  velocity  exactly  by  computing  the  Fourier  coefficients- 
of  the  data  over  the  whole  mission  (assuming  the  orbit  to  be  nearly  period¬ 
ical),  and  using  these  coefficients  to  calculate  the  acceleration,  since 
the  acceleration  coefficients  are  related  to  the  velocity  ones  by  the 
simple  relationship 


((accel.)j 


j(veloc) . 


The  Np  -  vector  of  acceleration  values  a  can  be  obtained,  formally, 
by  multiplying  the  Np  -  vector  of  velocity  values  by  a  NpxNp  "dif¬ 
ferentiator"  matrix  5  : 


a  =  S  vi; 


(4-1) 


The  model  for  the  accelerations  is,  therefore, 
a  =  S  A  c 


(4.2) 


and  the  "observed"  accelerations,  which  consist  of  differentiated  signal 
and  noise,  are 

-(obs)  =  S  -12(obs)  =  S  A  c  +  S  n 

The  least  squares  adjustment  estimator  is,  therefore  (paragraph  (2.8)) 
c  =  (AT  ST(0 '  )_1S  ApAT  ST  (0 '  )-1S  Vxz(obs) 


(4.3) 


(4.4) 


where  O'  is  the  variance-covariance  matrix  of  the  "differentiated" 
noise 

O’  =  E  { (Sn ) (Sn ) T }  =  S  E  {nnT}  ST 


*  S  0  S' 


(4.5) 
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(4.6) 


Replacing  (4.5)  in  (4.4) 

£  >  (aWiV^aP  ast  (s*)^d"1s*s  y1!(obs) 

■  (AV'A)-  ad-  v12(obs) 

where  "St"  indicates  "The  pseudo  inverse  of  S  "  (StSA  =  A)  ,  so  the 
least  squares  adjustment  estimates  of  the  coefficients  based  on  the  dif¬ 
ferentiated  data  are  identical  to  those  obtained  from  the  data  directly. 

The  errors  in  the  estimates  being  the  same,  any  conclusions  derived  for 
accelerations  are  identical  to  the  corresponding  conclusions  for  velocities. 
The  same  principle  applies  to  least  squares  collocation,  as  it  can  be 
seen  by  replacing  (A>ST(0 ' )-1S A)  with  (A  S’(D ‘ )“ 1S A  +  C-1)  in  the  reasoning. 

Figure  4.1  below  shows  the  true  position  of  the  satellites.  Si  and 
S2  ,  and  their  "pseudo-positions"  on  the  mean  sphere  S{  and  S2  (the 
"prima"  symbols  indicate  the  counterparts  of  real  elements  on  the  mean 
sphere).  The  relative  positions  shown  have  been  chosen  for  pure  convenience 
of  representation:  the  satellites  can  be  both  above,  both  below,  or 
on  either  side  of  the  surface  of  the  sphere. 


Figure  4.1:  True  Positions  in  Orbit  and 

Pseudo-Positions  on  the  Mean  Sphere 
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-  ***  - 


Let  / 

I  Co  ar{r *  .<*> '  .A  * )  -  r0  ar(r,<M)  j 

5a(Q‘,Q)  =  a(Q')  -  a (Q)  =  |  <£d  (r ‘  ,<p  1  ,A  ' )  -  <J>0  a^(r,4>,A)  (  (4.7) 

(id  a^(r'  ,<J> '  ,A ' )  -  A0  ax(r,$tA)J 

be  the  variation  in  inertial  acceleration  between  points  Q  and  Q' 
where  the  r  and  <p  components  have  been  modified  by  the  removal  of 
a  contribution  from  the  zonals  (paragraph  1.3).  As  it  happens,  this 

contribution  comes  only  into  ar  and  a  ,  a.  is  free  from  this  zero 

frequency  problem.  Let  *  A 

e  =  e( 2  ■  § 1 2 

be  the  difference  between  the  unit  vector  in  the  direction  of  the  actual 
line  of  sight,  e12  .  and  the  unit  vector  e{2  for  the  South-North  per¬ 
pendicular  to  the  radial  direction  at  the  "pseudo-midpoint"  P‘  .  The 
difference  between  the  true  line  of  sight  value  and  the  corresponding 
"pseudo-observation"  on  the  sphere  is  then 

6ai2  =  a12  -  &i2 

=  si2(a(Si)  -  a(S2))  -  ei2(a(S!)  -  a { S 2 ) ) 

=  ell  <5a(S(,S2)  -  (§12  -  <$a(Si,  S2) 

■  eiI(5a(S{,  Si  )  -  6a(Si,  S2))  +  e  6a(Si,  S2)  (4.8) 

Let  dmax  be  the  maximum  possible  distance  between  a  satellite  and  the 
corresponding  "pseudo-satellite"  position  on  the  sphere.  The  value  of 
dmax  will  depend  on  the  radial  and  horizontal  components  of  the  separation 
vector,  but  one  could  well  argue  that  a  purely  radial  separation  should 
have  much  the  same  effect  as  the  total  displacement  if  both  have  the 
same  size,  because  the  lateral  fluctuations  are  nearly  one  order  of  magni¬ 
tude  smaller  than  the  radial  one  (appendix  A).  Consider  the  first  term 
of  (4.8).  The  contribution  of  the  nth  harmonic  to  this  term  is 

61I  ( ( S { ,  S*)  -  6an  (Si,  S2)) 

where 

arn(Si)  rJ0-  arn(S2)  rjz)  Sa^  (Sa,  S2) 

6an(Si,  S2)  =  a^Si)  $0^-  a^n(S2)  $0  ^  =  <Sa^n  (Si,  S2)  (4.9) 

‘x„<s>>  J*'’  5*Xn  <S-  S>> 

arn  ,  a.n  ,  axg  being  the  sum  of  all  terms  of  degree  n  in  the  expansions 
of  ar  ,  Sa  ,  ax  (expressions  (1.8,  a-c)).  The  variation  of  arn  with 
radial  distance  is 

arn(R><j>»X)  -  arn(R  +  b,<+i,\)  s  (1  -  (^_)n+2)  arn(R,$,A)  (4.10) 
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and  similarly  for  a^n  ,  and  axn  .  As  both  satellites  follow  the  same 
orbit  with  a  separation  that  is  small  compared  to  the  wavr>ngth  of  the 
second  zonal  and  its  disturbances,  one  can  regard  them,  i  first  approx¬ 
imation,  as  moving  up  and  down  simultaneously,  so  their  h^.ghts  are  the 
same.  In  such  case 

'Sarn(Si ,  S2 )  -  <5arn(Si  ,  S2 )  s  (1  -  2)  6arn(S{  ,  S2) 

and  similarly  for  the  other  components  of  an  ,  so 

Sin(S{  ,  S2 )  "  ^3n(Si  ,  S2)  s  (1  -  )  ^?n(Si  >  S2 )  (4.11) 


The  magnitude  of  the  first  term  of  (4.8)  is,  therefore, 

eJ(6an(Si,  S^)  -  «in(Sx ,  S2))  <  1 1  e{2  j  |  (1  -  (^)n+2 )  M«§n(Si,  Si)\ 


As 


I  ei2|  |  ■  1  and  |  1 6an(Sx; ,  S2  )|  |  <  |j6an||  (where 


the  maximum  size  of  the  magnitude  of  an  oh  t¥lixmean  sphere),  and’Ifie  seg¬ 
ment  5’  has  constant  length,  it  follows  that 


|6aJlmtfie 


is 


ej(6 S2 )  -  5|n(Si ,  S2))  s  (1-  (^r4)  I  l«in(Si.  Si ) ]  ]  ma/  (4.12) 

The  difference  between  6a(Sx,  S2)  and  5a(Si,S2)  is  not  very  large,  so 
the  second  term  in  (4.8)  Ts,  approximately, 

eT  6a(Sx ,  S2)  s  gT  6a(Si ,  S2 ) 


R  ,n+  2, 


The  maximum  value  is 

J  6|(S i.  Si)  s  II § N max 


1 

'max 


The  magnitude  of  ||e||m  is,  to  a  first  approximation  (disregarding 
the  curvature  of  the  Eartn),  function  of  the  relative  displacement  of  the 
sate! 1 ites: 

||e||  s  ZSII  t  (4.13) 

1 1 -1 ‘max  pi2 


where  Ar12  is  radial,  Api2  along  the  line  of  sight,  and  Ati2 
across-track.  According  to  Appendix  A  ,  the  largest  displacements  due  to  the  an- 
amalous  field  are:  Api2max  =  0.7  km,  Ar12max  =  0.7  km,  Aii2!nax=  o.4km.  Ar  = 
3.9  km  (vertical  displacement  of  each  satellite).  Adding  1  km  to  Api2[Tia™x 
Ari2max  and  A'izmax  ,  to  introduce  an  extra  error  margin,  the  value 
of  iTsllmax  according  to  (4.13)  is 


The  contribution  of  the  nth  harmonic  to  the  second  term  of  (4.8),  therefore, 
is 


6an(S{ >  S2)  < 


0.009  lisa. 


'max 


(4.14) 
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(4.15) 


From  (4.8),  (4.12),  and  (4.14)  one  gets,  finally, 
^<-  [(1-  (■&?*)  +  0.009  n«in!lmax 


as  an  upper  bound  for  6ai2  •  The  maximum  value  of  h  is  d^^  ,  as 
explained  earlier.  This  maximum  absolute  displacement  of  any  or  the 
two  satellites  is 


(4.16) 


iirTOx  •  AT‘=«iax 

in  [n.  lb;  with  tneir  values  according  to  appendix  A  : 


(1  km  added  to  each) 


d  s  4.7  km 
max 


and  the  upper  bound  for  fa  12  becomes 

\n+2i 


fa12n  1  [1-009  -  (0. 99928087 )n+2]  li^nllmax 


"  pn  1 1 I  max 


(4.17) 


This  expression  limits  the  difference  between  the  actual  line  of  sight 
acceleration  along  the  non-circular  orbit,  and  the  corresponding  pseudo¬ 
observation  on  the  sphere.  This  is  the  largest  possible  value  (if  the 
argument  behind  (4.17)  can  be  accepted)  of  the  correction  that  has  to 
be  introduced  in  order  to  transform  the  original  data  set  into  the  corres¬ 
ponding  set  of  "accelerations"  on  the  mean  sphere.  The  factor  pn  is 

pn  =  [1.009  +  (0.99928087) 0+2 ]  (4.18) 


If  no  correction  is  applied  to  the  accelerations,  and  they  are  analyzed 
as  if  they  were  already  on  the  mean  sphere  and  on  the  ideal  orbit,  the 
error  in  the  estimated  coefficients  will  have  two  components:  (1)  that 
due  to  the  data  error,  which  has  been  listed  in  the  previous  section; 

(2)  the  error  due  to  the  fact  that  the  uncorrected  accelerations  are 
not  on  the  ideal  orbit.  As  both  errors  have  quite  independent  sources, 
the  total  rms  error  per  degree  is 

aen  =  Coen(M  +  aGn(2)]i  (4*19) 


where  cen(i\  corresponds  to  the  first,  and  ce  /2i  to  the  second  source 
of  error.  Tne  relative  error  per  coefficient  is,  then 


_  acn(2n+l)*1 
°n(T)  "  o7X2n+i w 


(4.20) 


where  Pn(n  is  the  propagated  noise  listed  in  Tables  3.1  or  3.2  of  section 
3.  If  the  only  difference  between  the  real  and  the  ideal  orbits  were 
an  increase  in  radial  distance,  orexpansion,  so  R  -  R'  ■  d^  were  constant, 
then  it  would  follow  from  (4.8),  (4.10),  and  (4.17)  that 
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(4.21) 


oai^ 


=  pn  ai2n 


as  ei2  =  ii2  and  6ax2n(Si,  S2)  - 
points  along  the  orbit,  so 


V) 


|5a12 


,  R 


max 


)  6a12n(S{,  Si)  at  all 


n1  max 


I  _1 
'  max 


=  Pr 


(4.22) 


In  reality  6ai2n  becomes  modulated  by  the  changes  in  distance  between 
true  and  ideal  orbit,  as  this  distance  cannot  be  constant.  Consequently 
expression  (4.21)  is  no  longer  true,  but  merely  an  approximation  to  a 
very  complicated  relationship.  The  error  should  be  less  than  with  a 
constant  radial  increase  dmax  at  almost  every  point,  so  perhaps  the 
use  of  the  right  hand  side  of  (4.22)  as  an  upper  bound  for  cn/.\ 
a  realistic  choice:  on(2)  s  Pn  • 


is 


coefficients  estimated  from  the  original  set  of  accelerations 
lied  Cftm(°)  ,  and  let  a12(0  =  al2(°)  +  6a12(°)  be  the 

(°)  of  6a j 2  computed 


Let  the 
ai2(°)  be  ca 

original  accelerations  correc^e 


ed  with  values  5ai2’ 

be  the  coefficients  obtained  from  the 


using  the  C .Let  CfLO  . 

analysis  of  the  ai2i  /  ,  and  let  6ai2l  >  be  a  new  estimate  of  £a1?  ,  . 

based  on  the  C$m  O'  .  Correcting  the  original  ax2(3)^  with  the^  5ai01/ 
and  then  analyzing  again,  one  obtains  .  If  /va12(2'  =  a12(a)+£aj3)  i; 

closer  to  the  transformed  set  of  pseudo-obs  ai2f"' 


subject  of  this  paragraph!  than  the 
be  an  improvement  pn  Cnm^1'  .  6a12(2) 

Let  C$m(M)  ,  6aiM',  P n  '  ( 2 ) ,  etc.  be  tl 


( 0 1  +  6a ; 2 
()/.,  then 
(0 


a  j2 
on  6a  i2 

«-*- >.  ^nm *  '  *  x c  ,  ci.\,.  uc  the  outcome  of 

analysis^correction-ana lysis  p^edure.  _As  M  -*•_ »  £,a 


(that 

Lnm 

o 


a] 


s  the 


nm 


nm;  —vi 

so  onvij*eon(i)' 
iterations  of  this 


may  conve 


to  some  limit  and  Pn(2)iM; pff(2) 

because  C9K  is  better  than  CnmvM).  As  Dn  is  the  ratio  of  6§1:n  toM  . , 
J3)  ,  if  the  relative  error  in  the  CSm(M-l)  ,  and  thus  in  the  6a12(m‘  { 


ai2n' 

is 

6ai2 


or 


so 


riiM-1)  .  the  error  ratio  for  the  corrected 


and  thus  for  the 

(«i 


where,,",Pn(2)  <  Dn 
ati  ^ 
thu 
(M) 


ai2 


he  6c. _ 

=  aj  2 ( 0 )  + 


{pn(M 


Pn  hi(T) 


(M-l)jJ 


s(M) 


■(H) 


*  P 


n(2) 


P2(P2 

rnvtn 


P2  (M-2) 
n(T) 


<  o2 

1  %  °n  (T) 


+  pn(1)) 


,(M-3) 


(4.23) 


qnpn(3)  ■  qn  P(T)  +  %  p2n(2)  +  pn  pn(x )  +  pn(x ) 


where  q 


n 

.a(H> 

n(T) 


n 

s  P, 


Repeating  this  reasoning  backwards 
(M-l) 

Ml) 


(pn(1) 


a  ]  +  Q  M-2)  o. 

Mn'  Hn 


1  times: 
•  +  q. 


n 


’n(l)  pn(1 ) 


M 


+  pf/i>  Cq, 


n(2) 


<M-D 


+  q, 


(M-2) 


+  qn  +  i] 


=  q«  +  p 


“07 


The  noise  Sn 

is  always  the  same,  so 


-  1-  9rP 

no)  rr*%i 

in  the  corrected  accelerations  ai2^  -  ajjp 


°n  ( 1 ) 
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is  independent  of  M 


(4.24) 


rtf!-” 


so 


co  Lim 
°n(T)  ~  M-*30  pn(T) 


if  Pn  <  1 


(M)  - 


Examining  (4.18)  one  can  see  that  for  n  =  331  (top  of  the  band),  is 


pn  «  0.22 


As  both  Pn  and  fall  with  decreasing  n  ,  for  all  n  in  the  band  is 

0n(T)(H>  •=  p”‘(T)<M>  and  pn  <  "in  <  1 
Therefore  it  is  true  that 

pn(T)  S  TT <4'2 

The  values  of  Pn(T)  and  of  On(l)  have  been  listed  below,  together  with 
the  corresponding  values  of  the  error  Pn(x)  given  in  Table  3.2  and 
derived  under  the  simplifying  assumptions.  ;The  closeness  between  the 
"true"  and  the  "ideal"  relative  errors  is  clear.  Moreover,  the  iterative 
procedure  appears  to  converge  very  quickly,  as  Pn(T)  and  pn(T)  4 
agree  to  three  decimal  places.  ^ 


n 

pn(l)  x  100 

pn(T)  *  100 

pn(T)  x  100 

10 

0.0063 

0.0063 

0.0063 

100 

0.0686 

0.0688 

0.0689 

150 

2.8401 

2.8582 

2.8583 

200 

7.9146 

7.9982 

7.9987 

250 

21.620 

21.9581 

21.9583 

300 

65.941 

67.3613 

67.3614 

331 

82.269 

84.3747 

84.3748 

It  must  be  remembered  that,  on  the  one  hand,  the  reasoning  leading  to 
(4.25)  is  by  no  means  rigorous,  while, on  the  other  hand,  the  assumptions 
made  in  it  have  been  mostly  on  the  conservative  side.  It  appears  that 
a  series  of  iterations  as  described  here  should  transform  the  original 
set  of  observations  into  a  set  of  "accelerations"  on  the  ideal  orbit, 
in  such  way  that  the  coefficients  estimated  from  these  transformed  set 
would  have  much  the  same  accuracies  as  those  listed  in  section  3  and 
derived  under  the  simplifying  assumptions  of  paragraph  2.1.  In  this 
sense,  the  results  of  section  3  are  supported  by  those  above;  they 
may  very  well  indicate  the  maximum  amount  of  information  on  the  geopotential 
that  can  be  extracted  by  linear  estimation  techniques  from  a  real  set 
of  SST  data. 


4.3  The  Effect  of  Errors  in  the  Calculated  Orbits 


To  apply  the  corrections  5axiM)  to  the  a, 2  ( 3 )  one  must  know 
the  position  of  the  two  satellites,  Sx  and  S2  .  In  the  previous 
paragraph  it  was  assumed  that  this  knowledge  was  exact .  In  reality  , 
however,  the  gravity  field  model  used  to  compute  the  orbits,  the  coordin¬ 
ates  assigned  to  the  tracking  stations,  and  the  tracking  d<->*a  used  to 
calculate  the  orbits,  all  contain  errors,  so  there  is  a  discrepancy  between 
computed  and  true  orbits.  With  present  tracking,  coordinates,  and  models, 
the  errors  in  the  calculated  orbits  are  not  likely  to  exceed  a  few  meters. 
The  direct  effect  of  these  errors  is  an  additional  error  in  the  corrections 
5a12(M)  ,  as  the  accelerations  aJ2(°)  +  5a12(M)  are  "dropped"  from 
their  true  positions  to  the  calculated  ones.  The  change  in  cai2,y  ' 
for  a  displacement  of  As  meters  is  likely  to  have  a  size  comparable 
to  the  effect  of  a  purely  vertical  displacement  of  the  same  magnitude. 

From  (4.11)  follows  that  such  change  is 


AAa  <M)  -  M  (  R 

A6a^n  (1  (7TTIs 


(4.26) 


or,  to  a  very  good  first  approximation 


Aa12 


(M) 

n 


=-(n  +  2) 

=  6n(As)  Sa12 


(M) 

da12n 


(M) 


(4.27) 


Even  *nr  n  =  331  ,  at  the  very  top  of  the  band,  the  relative  error 
3r  caused  by  an  orbital  error  of  10  m  is  insignificant:  S33i(a0'  = 

5x  10~ s.  This  error  decreases  with  n  ,  as  indicated  by  (4.27),  so  it 
is  smallest  at  n  =  0  .  However,  the  size  of  the  zero  harmonic  is  so 
large  that  its  net  change  over  10  m  amounts  to  some 
3  mgal  .  If  no  correction  term  for  this  and  the  other  even  zonal s  is 
introduced  in  the  definition  of  the  derivative  of  the  line  of  sight  velocity, 
as  suggested  in  section  1  ,  then  the  effect  of  a  vertical  displacement 
is 


So (As)  fa 


12  0 


(m)  _  9  .  \b_  3  GM 

0  - *2  Sin  -*R-  - — 7 

2  3  r  r 


AS 


(4.28) 


or  0.138  mgal  every  10  m  .  Because  the  orbital  error  is  not  constant, 
this  change  in  the  zero  harmonic  constribution  is  modulated  by  the  complex 
shape  of  the  error,  so  the  power  of  the  zero  harmonic  change  spills  over 
the  whole  spectrum  corrupting  all  the  estimated  coefficients.  As  the 
power  associated  with  n  >  100  ,  in  mgal,  is  well  below  0.138 
mgal  at  satellite  altitude,  this  can  have  serious  consequences.  In  the 
model  adopted  here  for  ai2  the  zero  harmonic  has  been  excluded,  so 
only  the  errors  due  to  the  other  terms  in  the  expansion,  all  of  them 
quite  negligible,  are  present. 
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The  objective  of  this  section  is  to  include  some  thoughts  on  how  to 
process  the  masses  of  information  collected  from  a  SST  experiment  into 
a  global  geopotential  model  of  very  fine  resolution.  This  problem  goes 
beyond  the  scope  of  this  study,  but  some  further  elaboration  of  concepts 
introduced  in  previous  sections  may  help  clarify  this  difficult  subject, 
and  perhaps  point  out  directions  for  future  research. 


5.1  An  Iterative  Approach 

The  argument  used  in  the  previous  section  to  substantiate  the  num¬ 
erical  results  of  section  3  suggests  that  a  model  may  be  obtained  through 
successive  approximations.  To  be  able  to  use  the  sort  of  "analysis-correc¬ 
tion-analysis"  approach  hinted  at  in  that  argument,  it  is  necessary  to 
have  certain  regularities  in  the  data  that,  though  physically  possible, 
may  not  occur  in  practice.  On  practical  grounds  one  can  question  two 
assumptions  made  implicitly  in  last  section:  that  the  deoartures 
from  perfectly  circular  orbits  were  due  to  the  anomalous  gravity  field  alone, 
and  that  the  data  were  sampled  uninterruptedly  at  constant  intervals  during 
the  whole  mission,  all  measurement  errors  having  the  same  standard  deviation. 

Even  when  the  compensating  mechanism  could  eliminate  all  non-gravi- 
tational  forces,  and  when  gravitational  fields  (other  than  the  Earth's) 
and  the  effect  of  the  body- tide  could  be  calculated  exactly  from  existing 
models  and  thus  discounted,  the  relative  positions  of  the  satellites  are 
not  going  to  vary  as  predicted  from  the  action  of  the  anomalous  field 
alone.  The  reason  for  this  is  that  it  is  impossible  to  determine  and 
control  exactly  the  state  of  each  satellite  at  "injection"  time,  to  make 
sure  that  both  move  along  the  same  orbit,  and  that  such  orbit  is  as  close 
to  circular  as  the  field  would  allow.  As  a  consequence, besides  the  relative 
motion  caused  by  the  field  there  will  be  a  "drift"  due  to  incorrect 
initial  conditions.  This  drift  will  follow  a  more  or  less  arbitrary  dir¬ 
ection,  resulting  in  the  spacecrafts  moving  towards  or  away  from  each 
other  until  their  separation  and  the  direction  of  the  line  of  sight  change 
so  much  that  observations  taken  at  different  times  cannot  be  described 
sufficiently  well  by  equations  that  assume  a  fixed  distance  and  angle.  Since 
the  total  relative  motion  can  be  calculated  from  the  reference  orbits 
with  an  accuracy  of  a  few  meters  ,  the  change  in  relative  configuration 
with  time  should  soon  become  apparent.  This  change  should  have  a  more 
or  less  periodic  component  due  to  the  irregular  gravitational  field-and 
a  trend  due  to  previous  orbital  manouvers'  errors.  While  these  errors 
cannot  be  determined  accurately,  the  drift  itself  increases  with  time 
until  it  becomes  easy  to  estimate.  The  fact  that  the  satellites  are 
tracking  each  other  continuously,  in  addition  to  being  tracked  from  ter¬ 
restrial  stations,  should  help  to  obtain  a  good  estimate  of  the  drift 
and  the  drift-rate.  When  the  drift  reaches  a  maximum  allowed  value,  the 
controlling  rockets  of  either  spacecraft  can  be  fired  briefly  to  reverse 
the  drifting  motion,  so  the  drift  begins  to  decrease  (the  reversal  need 
not  be  exact).  This  scheme  is  a  "dead  zone"  control  policy  where  correcting 
action  is  applied  only  when  a  certain  limit  is  about  to  be  exceeded.  As 
the  purpose  is  to  reverse  the  drift  rate  rather  than  to  eliminate  the 
drift  altogether,  the  correction  need  not  be  particulary  drastic.  Even 
so,  the  use  of  such  a  scheme  must  increase  the  amount  of  fuel  required 


during  the  mission,  and  thus  the  overal1  cost.  Of  course,  some  cor¬ 
recting  maneuvers  of  the  sort  described  here  will  have  to  take  place  from 
time  to  time,  to  stop  the  two  spacecraft  from  separating  too  much  (so 
the  radar  beam  does  not  cut  too  deeply  inside  the  atmosphere,  or  is  inter¬ 
cepted  by  the  curvature  of  the  Earth  itself,  at  a  separation  of  about  2700 
km)  or  from  becoming  too  close,  as  the  shorter  their  distance,  the  weaker 
the  signal  and  the  less  accurate  the  results,  as  shown  in  section  3,  Tables 
3.6,  3.7,  and  3.8.  The  question  is  how  often  such  maneuvers  will  take 
place. not  whether  they  shall  be  needed  at  all.  Compensatory  changes  in 
the  design  of  the  mission  could  allow  the  use  of  "dead  zone"  control  to 
maintain  relative  orientation  without  raising  the  cost  coo  much.  It  may 
be  possible  to  choose  a  higher  orbit,  where  less  fuel  for  drag  compen¬ 
sation  is  needed  (the  decrease  in  fuel  requirement  with  height  should 
be  quite  fast)  so  the  balance  can  be  used  for  maneuvering.  The  length 
of  the  mission  could  be  increased,  the  accuracy  of  the  data  improved, 
or  both,  to  compensate  for  the  weaker  signal  at  greater  height.  All  of 
these  factors,  some  contradictory,  and  others. as  well,  can  only  be  bal¬ 
anced  properiy  in  a  thorough  mission  design  study  incorporating  the  demands 
of  data  processing  among  the  main  questions  to  be  considered.  Given  the 
magnitude  of  such  demands.  SGme  regard  for  them  appears  natural. 


Another  assumption  made  in  section  4  that  may  not  be  fulfilled 
in  practice  is  the  existence  of  uninterrupted  series  of  measurements  lasting 
the  whole  mission.  Breaks  will  occur  in  the  data,  some  much  too  large 
to  be  ignored  or  "patched  up"  by  interpolation  from  surrounding  data. 

There  may  be  fluctuations,  as  well,  inthequality  of  those  measurements 
due  to  problems  in  the  satellites  themselves,  or  in  the  surrounding  medium, 
as  in  the  case  of  severe  ionospheric  perturbations .  So  the  stream  of 
data  will  not  by  perfectly  uniform  in  quality  (variations  in  the  standard 
deviation  of  the  errors)  or  unbroken.  These  departures  from  the  assump¬ 
tions,  if  severe  enough,  would  make  a  close  application  of  the  ideas  presented 
in  earlier  sections  quite  impossible.  One  could  begin,  as  an  alternative, 
by  differentiating  numerically  the  data  to  obtain  line  of  sight  accelerations. 
Supposing  that  the  model  of  section  1  is  correct.,  these  accelerations 
on  a  sphere  have  the  form 
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where,  according  to  expression  (2.9)  and  (2.15,a-b) 
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Functions  that  can  be  expanded  in  series  of  the  type  of  (5.1),  which 
includes  the  ordinary  spherical  harmonic  expansion,  can  be  analyzed  to 
obtain  the  values  of  the  by  very  efficient  methods  resembling 

numerical  quadratures.  The  derivation  and  implementation  of  optimal  methods 
that  minimize  the  estimation  error  in  the  presence  of  unhomogeneous  and 
and  correlated  noise  have  been  discussed  by  Colombo  (1981,  paragraphs  2.9 
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and  2.12).  While  the  procedures  presented  by  this  author  in  that  work  have 
been  derived  for  the  case 


(i.e.,  spherical  harmonics)  the  extension  to  the  more  general  case  of 
(5.1)  is  quite  direct.  The  only  requirement  is  that  the  data  occupy  the 
nodes  of  a  regular  grid  where  the  separation  between  meridians  is  constant. 
This  can  be  done  by  interpolating  the  differentiated  velocities  on  those 
nodes  from  surrounding  daca  points,  while  also  determining  che  correlations 
and  standard  deviations  of  the  interpolated  values  fromthose  of  the  obser¬ 
vations,  as  these  ure  needed  to  set  up  the  optimal  estimator.  If  no  reliable 
interpolation  is  possible  on  some  node,  because  of  large  data  breaks  in 
the  vicinity,  a  value  of  zero  with  a  “standard  deviation"  ecual  to  the 
rms  of  the  line  of  sight  acceleration  could  be  used  instead. 

By  setting  up  an  appropriate  control  mechanism  to  maintain  the  relative 
conf iguration  of  the  satellites  within  sufficiently  close  limits,  and 
by  using  the  interpolated  accelerations  as  pointed  out  above,  the  interative 
scheme  could  proceed  basically  along  the  lines  of  section  4,  except  that 
the  data  would  be  analyzed  as  in  Colombo  (1981),  rather  than  as  in  section 
2.  If  such  a  control  scheme  proves  to  be  feasible,  the  next  important 
question  i§  how  to  calculate  the  corrections  6&n  of  expression  (4.8) 
from  the  C'nm  estimated  in  the  previous  iteration,  to  refine  the  pseudo¬ 
observations  on  the  mean  sphere.  Because  of  the  irregular  shape  of  the 
orbit,  calculation  of  these  corrections  using  exact  relationships  is 
too  laborious  to  be  practical,  even  with  very  powerful  computers.^  The 
main  reason  is  the  rumber  of  operations  needed  to  obtain  every  can  > 
which  is  proportional  to  the  number  of  coefficients  (some  11  x  10“  ** 

N  =  331).  This  matter  needs  thorough  investigation,  but  the  answer  should 
involve  some  sort  of  approximation,  to  reduce  the  computer  burden.  One 
possibility  that  may  be  worth  exploring  is  as  follows:  consider  a  spherical 
shell  extending  5  km  above  and  5  km  below  the  mean  orbital  sphere.  The 
whole  orbit,  according  to  Appendix  A,  should  lie  within  this  shell.  Imagine 
the  shell  subdivided  according  to  a  regular  grid,  equal  angular  for  instance, 
where  the  blocks  are  yD  x  y°  in  size,  so  the  shell  is  partitioned 
into  cells  each  yJ  x  y°  x  10  km  in  volume  (the  size  of  y  should  be 
decided  by  detailed  study).  The  vertices  of  the  cells  are  arranged  in 
equal  angular  fashion,  so  one  can  compute  the  three  inertial  acceleration 
components  at  every  vertex  according  to  expressions  (1.8,  a-c)  very  effic¬ 
iently  using,  for  example,  algorithms  like  those  described  by  Colombo 
(1981).  To  calculate  the  correction  5ax2  one  must  know  the  value  of 

in  the  actual  orbit  and  in  the  ideal  orbit  ,  and  for  that  the  three 
accelerations  at  points  Si  and  S;  ,  S{  and  S)  (in  Fig.  4.1)  are 
needed.  These  accelerations  can  be  computed  from  their  values  at  the 
vertices  of  the  cell  containing  those  points  by  some  interpolator  pro¬ 
cedure,  which  3nould  be  a  great  deal  easier  than  an  exact  calculation. 

The  type  of  approximate  calculation  of  the  acceleration  components 
just  described  may  be  used,  as  well,  to  obtain  the  forcing  function  due 
to  gravitation  needed  to  integrate  numerically  satellite  orbits  with  a 
very  high  degree  and  order  field'  model,  like  the  one  whose  determination 
is  being  discussed  here.  Such  calculations  may  be  necessary,  for  example, 
to  reduce  the  orbital  errors  at  each  step  of  the  iterative  modelling  pro¬ 
cedure. 
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5.2  Other  Methods 


Instead  of  converting  the  data  into  Dseudoobservations  on  the  mean 
sphere  by  successive  iterations,  this  could  be  done  directly  by  estimating 
the  pseudoobs.  in  one  step  from  neighbouring  data  points.  This  means 
carrying  out  a  series  of  local  reductions  that  finally  covers  the  whole 
sphere  with  a  regular  grid  of  estimated  values.  These  values  should  cor¬ 
respond  to  a  variable  that  can  be  described  by  an  expansion  of  the  type 
of  (5.1),  and  which  does  not  have  to  be  the  line  of  sight  acceleration 
of  equation  (2.24).  These  local  reductions  can  be  done  \n  many  different 
ways.  One  possibility  is  to  use  least  squares  collocation  in  a  local 
manner.  A  problem  with  local  solutions  by  collocation  is  thatthe  signal 
in  SST  measurements  tends  to  be  too  strongly  correlated  over  considerable 
distances,  and  this  causes  the  variance-covariance  matrix  to  be  too  ill- 
conditioned  to  be  inverted  in  a  dependable  way.  According  to  R.  Rummel 
(private  communication)  a  minimum  separation  of  about  100  km  between  data 
points  is  necessary  for  stable  inversion.  As  the  actual  points  are  likely 
to  be  spaced  some  30  km  apart  along-track  (with  a  sampling  interval  of 
4  s  )  and  by  less  than  that  across-track,  (because  the  separation  between 
adjacent  passes  should  be  of  only  a  few  kilometers  by  the  end  of  a  six- 
month's  mission)  some  kind  of  paring,  or  decimation,  of  the  data  will 
be  needed  to  achieve  the  larger  spacing.  This  could  be  done  in  the  first 
of  two  steps:  beginning  with  decimated  data,  a  set  of  pseudoobservations 
could  be  obtained  without  stability  problems,  and  a  first  estimate  of 
the  harmonic  coefficients  could  be  made  from  this  set.  In  the  second 
step,  all  the  data  could  be  used  after  substracting  fromthem  their  nominal 
values  according  to  the  model  produced  in  step  one,  which  could  have  high 
degree  terms  already.  Such  residuals  are  going  to  be  less  correlated 
than  the  original  measurents,  because  much  of  their  low  frequency  content 
would  have  been  removed,  so  the  inversion  of  their  variance-covariance 
matrix  may  be  stable  in  spite  of  the  close  spacing  of  the  data  points. 

If  no  control  of  the  relative  alignment  of  the  satellites  takes  place 
and  they  are  allowed  to  drift  freely  with  respect  to  each  other,  the  cor¬ 
responding  observation  equations  would  lack  the  regular  structure  assumed 
in  section  2.  In  such  a  case  the  ideas  for  analysing  the  data  discussed 
in  the  previous  paragraph  are  not  applicable,  and  only  non-iterative  methods 
like  the  one  just  outlined  seem  to  offer  any  real  hope. 

Any  method  that  first  creates  a  regularly  spaced  set  of  pseudoob¬ 
servations  on  the  mean  sphere  and  then  analyses  it  efficiently  to  obtain 
the  potential  coefficients  requires  two  main  things: 

(a)  a  model  for  the  data  that  is  both  practical  and  accurate.  An 
example  of  this  may  be  the  model  adopted  in  sections  1  and  2,  but  this 
idea  needs  further  validation,  as  pointed  out  at  the  end  of  paragraph 
1.3; 

(b)  an  efficient  and  accurate  way  of  reducing  the  measurements  in 
the  actual  orbit  to  pseudoobservaticrs  on  the  mean  sphere. 

Both  problems  require  more  research  to  clarify  them,  and  this  clarifi¬ 
cation  may  be  absolutely  necessary  before  developing  an  effective  technique 
for  processing  SST  data. 
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5.3  The  Use  of  Local  Solutions 


As  mentioned  in  the  introduction,  there  have  been  studies  of  the 
SST  problem  that  have  concentrated  on  purely  local  solutions  that  use 
data  taken  from  inside  a  relatively  small  neighborhood  of  the  points 
at  the  Earth's  surface  where  certain  quantities  associated  with  the  grav¬ 
itational  field  are  estimated.  The  advantage  of  the  global  approach  is 
that  every  value  that  is  estimated  is  obtained  from  the  processing  of 
all  the  data  available,  so  it  can  be  more  accurate  than  a  local  solution 
based  only  on  part  of  the  data.  Global  solutions  may  be  free  from 
the  numerical  instabilities  that  tend  to  be  associated  with  local  solutions 
Global  solutions  have  also  some  important  limitations.  A  spherical  harmon¬ 
ics  model,  for  instance,  must  be  truncated  at  some  finite  degree  for 
practical  reasons  and  this  limits  the  size  of  the  finest  detail  that 
the  model  can  represent.  Even  if  the  number  of  coefficients  is  not  a 
problem,  some  very  fine  but  also  quite  strong  features  that  may  be  sensed 
by  the  satellite  pair,  such  as  anomalies  along  ocean  ridges  and  trenches, 
mountain  ranges,  etc.,  are  likely  to  be  smoothed  out  by  an  optimal  global 
estimation  procedure  because,  on  a  global  scale,  they  are  similar  to 
measurement  noise.  Such  small  but  marked  features  may  be  recovered  by 
using  local  estimation  methods,  applied  in  the  knowledge  that  sharp  field 
variations  may  occur  in  certain  areas.  These  methods  may  use  residual 
SST  data,  with  respect  to  the  global  solution,  to  ensure  numerical  sta¬ 
bility  and  the  removal  of  trends  of  non-local  nature. 

Local  modelling  should  be  regarded,  therefore,  as  complementary  to 
global  analysis,  because  its  careful  application  in  selected  areas  may 
push  the  level  of  resolution  to  the  very  limits  allowed  by  the  information 
contained  in  SST  data.  Local  solutions  may  permit  also  the  combination 
of  SST  data  with  terrestrial  measurements  of  gravity,  with  satellite  altim¬ 
etry  over  the  oceans,  and  with  knowledgeof  the  geological  structures 
that  may  be  responsible  for  some  of  the  high  frequency  content  of  the 
signal,  all  of  which  may  be  difficult  to  incorporate  into  a  olobal  solution 


6.  Conclusions 


According  to  the  theory  given  in  sections  1  and  2  ,  assuming  that 
the  power  spectrum  of  the  geopotential  is  as  described  in  paragraph  3.1, 
and  that  the  reasoning  of  section  4  is  valid,  the  results  presented  in 
section  3  can  be  summarized  as  follows:  with  two  satellites  in  nearly 
the  same  circular,  polar  orbit,  at  a  height  of  160  km,  300  km  apart,  with 
a  tracking  accuracy  of  /?  x  10-6  m  s*1,  a  sampling  period  of  4  s  ,  an 
averaging  period  of  4  s  ,  and  a  mission  length  of  six  months,  the  coef¬ 
ficient  s of  the  spherical  harmonic  expansion  of  the  potential  up  to 
degree  n  =  331,  may  be  estimated  to  the  following  entent: 

(1)  the  relative  accuracy  of  the  potential  coefficients  could  be 
better  than  1%  for  r  $  130  ,  than  10%  for  n  t  210  ,  and  than  50%  for 
n  s  270  ,  using  least  squares  collocation.  With  least  squares  adjustment, 
the  results  may  be  the  same  up  to  n  =  200  ,  and  for  higher  degrees 
col  location  may  work  better. 
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(2)  the  accuracy  of  the  geoid  undulation  implied  by  the  coefficients 
could  t>e  better  than  0.05  mm  rms  for  wavelengths  of  between  3000  km  and 
40030  km,  and  better  than  10  cm  rms  in  the  band  between  140  km  and  3000  km. 

The  method  of  analysis  of  section  2  can  be  used  for  studying 
missions  where  the  orbital  plane  is  oblique  to  the  equator.  This  would 
require  some  changes  in  the  programs  listed  in  Appendix  B  ,  subroutine 
0NEREV  in  particular,  because  those  are  written  for  the  case  of  polar 
orbits  only.  As  they  are,  they  can  be  used  to  carry  out  an  analysis 
of  the  magnitude  reported  here  in  any  modern  computer  w'th  some  1.5  megabytes 
of  core.  If  the  highest  degree  studied  is  N  =  331  ,  the  central  processor 
unit  time  required  is  of  the  order  of  one  hour.  The  effect  on  the  results 
of  changes  in  some  of  the  mission  parameters  can  be  studied,  in  the  case 
of  least  squares  adjustment,  using  the  results  of  section  3  and  the  simple 
formulas  of  paragraph  2.13.  Only  potential  coefficients  and  geoid  undul¬ 
ations'  accuracies  have  been  calculated;  the  accuracies  of  other  quantities, 
gravity  anomalies,  for  example,  can  be  obtained  in  a  simple  way  from  the 
coefficients'  accuracies,  a  complete  listingof  which,  decree  by  dearee, 
appears  in  Appendix  C  . 

The  method  of  section  3  allows  for  a  spherical ,  rotating  Earth,  and 
for  data  consisting  of  velocity  averages.  A  number  of  simplifying  assump¬ 
tions  listed  in  paragraph  2.1  are  substantiated  in  sections  1  and  4,  in 
an  attempt  to  show  that  the  results  listed  in  section  3  represent  the 
limit  of  accuracy  for  a  global  model  obtainable  from  real  SST  data  if 
the  mission  could  be  carried  out  without  a  single  fault. 

The  accuracies  listed  being  global,  the  undulation  errors  are  likely 
to  be  worse  in  some  areas,  perhaps  where  the  field  is  strongly  anomalous, 
such  as  mid-ocean  ridges,  ocean  trenches,  mountain  ranges,  etc.,  and  nec¬ 
essarily  better  than  average  over  the  remainder  of  the  planet.  In  those 
areas  where  the  model  may  perform  poorly,  local  solutions (using  the  model 
as  a  reference  field  to  calculate  the  residuals  of  satellite  and  terrestrial 
data,  and  these  residuals,  in  turn,  as  the  observed  values)  may  provide  the 
finer  detail  that  a  global  technique  alone  cannot  reveal. 

As  argued  in  section  5  ,  an  iterative  solution  based  on  the  ideas 
of  sections  2  and4may  be  used  for  the  actual  analysis  of  SST  data.  Alter¬ 
natively,  non-iterative  methods  could  be  developed  for  that  purpose.  All 
such  methods  should  have  in  common  the  fact  that  they  rely  on  some  conven¬ 
ient  aproximation,  as  rigorous  solutions  such  as  those  based  on  the  "numerical" 
or  the  "analytical "  (i.e.,  celestial  mechanics)  approaches  are  virtually 
impossible  to  implement,  given  the  enormous  number  of  harmonic  coefficients 
to  be  adjusted  and  of  observations.  Two  essential  tasks  to  be  accomplished 
through  further  research  before  a  satisfactory  technique  can  be  found 
are,  in  all  likelihood,  the  following: 

(a)  validation  of  the  model  of  the  line  of  sight  relative  acceleration 
used  in  sections  1  and  2  ,  or  its  replacement  by  another  that  provides 

a  better  approximation  and  is  just  as  tractable  mathematically; 

(b)  development  of  a  satisfactory  procedure  for  reducing  the  SST 
measurements  to  a  set  of  pseudoobservations  regularly  distributed 

over  the  mean  orbital  sphere,  whose  analysis  can  then  be  carried  out  with 
efficient  algorithms. 


Further  work  on  how  to  model  the  SSI  signal  should  also  help  to 
clarify  the  question  of  the  influence  of  orbital  errors  on  the  estimated 
quantities,  influence  that..,  according  to  the  model  adopted  here,  may  be 
small . 
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Appendix  A:  Orbital  Perturbations 

A  spacecraft  orbiting  a  planet  cannot  follow  a  perfectly  circular 
orbit,  as  assumed  in  section  2,  but  must  move  in  a  more  irregular  course* 
because  of  anomalies,  or  fluctuations;  in  the  actual  gravity  field  compared 
to  that  of  a  central  point  mass.  In  the  case  of  the  Earth,  most  of  this 
effect  is  due  to  the  attraction  of  the  equatorial  bulge.  In  the  spherical 
harmonic  expansion  of  the  anomalous  field,  the  influence  of  this  bulge 
is  represented  mostly  by  the  second  zonal.  This  term  is  much  larger  than 
all  the  others  in  the  expansion,  so  most  of  the  orbital  perturbation  is 
due  to  this  term  alone.  The  following  calculations  will  take  into  account 
the  second  zonal  in  detail,  and  for  the  rest  a di splacement  of  200  m  in 
every  direction  will  probably  account  well  enough,  provided  no  strong 
resonances  occur,  as  assumed  in  paragraphs  2.1  and  2.6.  For  relative 
displacements  between  satellites,  400  m  will  be  added  to  those  caused 
by  the  second  zonal  alone.  Another  assumption  made  here  is  that  the  force 
corresponding  to  the  second  zonal  is  the  same  along  the  true  orbit  as 
along  the  ideal,  mean  circular  orbit.  As  the  total  displacement  of  such 
craft  is  mostly  vertical  and  of  less  than  5  km,  it  is  easy  to  determine 
that  the  force  of  the  second  zonal  cannot  change  by  more  than  0.3%,  so 
it  may  be  acceptable  to  regard  it  as  having  the  same  value  along  either 
orbit.  As  the  rotation  of  the  Earth  has  no  effect  on  the  force  of  a  purely 
zonal  field,  it  will  be  ignored.  The  motion  in  the  orbital  plane  which 
is  perpendicular  to  the  equator,  is  supposed  to  be  periodical,  which 
is  mathematically  possible  with  the  right  initial  conditions. 

Consider  the  system  of  inertial  coordinates  with  axes  x  and  y  , 
where  f  coincides  with  the  polar  figure  axis  of  the  Earth,  and  1  is 
in  the  equator  (figure  A.  I).  The  polar  coordinates  r  and  O'  cor,  »spond 
to  a  point  moving  along  the  orbit  with  an  approximately  uniform  circular 


motion  represented  by 

=  ^MODULE  2tt 

The  relationships  between  the  components  of  the  inertial  acceleration 
in  the  (x  ,  y)  and  the  (r  ,  o')  systems  are 

ax  =  ar  cosO’  -  a^  sino'  (A.l.a) 

a„  =  a„  sino'  +  a,  coso'  (A. 1 ,b ) 

y  r  0 

or  ax(t)  =  cos,J‘'t  ‘  sinwt  (A. 2, a) 

ay(t)  =  ar(t)  sinwt  +  a^t)  coswt  (A.2,b) 

From  the  definition  of  the  unnormalized  Legendre  function  introduced  in 
paragraph  2.2,  the  polar  components  of  the  acceleration  are 

ar(t)  =  "3  fr  (T):  Lnm(wt>  (A.3,a) 

=  3  “T"  3  2  (7-  cos  2ut  -  ^-) 
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v*> =  *  (f )!  fela,t) 

=  T  ‘^T  a2  s^n^u)t 

where  J2  =  C20  (unnormalized) 


Figure  A.l:  Acceleration  and  Position  in  the  (x,y)  and 
(rx  ({i1  =  at)  Coordinate  Systems. 


The  radial  perturbation  in  the  position  of  the  satellite  is 
Ar(t)  =  Ax(t)  coswt  +  Ay(t)sina)t 

=  [  dt '  av(t")  dt"  +  Ax0t  +  Axo]  coswt 
0  0  x 

^  t 1  • 

+  [  JT  dt  *  /o  ay  (t';  dt"  +  Ay0t  +  Ay0]  sinwt  (A. 4 

where  Ax,Ax0  etc.,  are  initial  conditions.  As.  Ar(t)  is  supposed  to 
be  periodical  with  the  right  Ax0  ,  Axo  ,  Ay  ,  A y0  ,  and  ignoring  a  constant 

form  Mint  morolv  <  hnn<jo'.  (-ho  *  i  ;n  nf  (hcs  mhil  hy  n  Tow  I-  i  1  <«nn  ( m  .  I  hr 

radial  displacement  is 

Ar(t)  -  -r  J;  COS  2wt  (A. 5 


according  to  (A.2,a-b),  (A.3,a-b),  and  (A. 4).  Furthermore  u>  =  / ,  so 

Ar(t)  =  02  cos2wt  (A. 6) 

and  using  the  values 
a  =  6371  km 
r  =  6531  km 
J2  =  1.1  x  1(T3 
(A. 6)  becomes 

Ar(t)  =  -3.7  cos2wt  [km]  (A. 7) 

So  the  amplitude  of  the  oscillation  due  to  J2  alone  is 

ArMAX  =  3-7  km  (A. 8) 

The  relative  radial  displacement  between  two  satellites,  assuming  that  the 
change  in  their  distance,  and  thus  in  ip  ,  can  be  ignored,  is 

Ar12(t)  =  -3.7  (cos2wt  -  cos2(u)t  +  \p)) 

=  -3.7  x  2sinp  sin(2wt  +  ip) 

=  0.33  sin(2wt  +  ij;)  (A. 9) 

for  an  intersatell ite  distance  p  =  300  km  ,  so 

Ar12MAx  =  °-3^  km  (A. 10) 


The  relative  motion  along  the  line  of  sight  can  be  estimated  by 
integrating  the  relative  velocity  v12  .  As  this  velocity  is  nearly  a 
sinewave  in  time,  of  frequency  2w  ,  the  expression  for  the  variation 
in  p  is 

Ap12(t)  =  /?  V^IE-)cos(2aJt  +  6)  (A. 11) 

where  B  is  some  phase  angle  of  no  consequence  here.  Replacing  v2 2 ( rms ) 
with  its  total  value  according  to  Table  1.1,  -md  adopting  oj  =  1.2  x  10*3 
rad  s*1  (as  r  s  6531  km), 

AOl2MAX  =  0,3  km  (A. 12) 

Because  of  the  force  is  due  to  a  zonal,  the  across-track  displacement  is 

Ati2(t)  =  0  for  all  t  .  (A. 13) 
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The  absolute  along-track  displacement  of  each  satellite  can  be 
obtained  from  an  expression  similar  to  (A. 7),  but  it  is  not  needed  in 
section  4.  Finally,  adding  200  m  or  400  m,  as  the  case  may  be,  to  account 
for  the  rest  of  the  anomalous  field,  the  perturbations  amount,  approximately 
to 


ArMAX  =3. 7+0. 2  =3. 9  km 
Ari2j^y  =  0.3  +  (2  x  0.2)  =  0.7  km 
ADl2MAX  =  +  (2  x  0,2)  =  0.7  km 
At12max  s  0.0  +  (2  x  0.2)  =  0.4  km 
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Appendix  B:  Computer  Programs 


This  Appendix  contains  descriptions  and  listings  of  the  main  programs 
and  subroutines  used  for  the  error  analysis  whose  theory  is  contained 
in  sections  1  and  2  and  the  results  of  which  appear  in  section  3. 


B.l  Main  Program 

Two  main  program  versions  were  developed:  the  first  one  creates 
the  non-zero  diagonal  blocks  of  (AT  D_1  A)  ,  stores  them  in  unit  10  (a 
tape  file),  creates  the  normal  matrix's  blocks  by  adding  C-1  or  not, 
depending  on  whether  least  squares  collocation  or  least  squares  adjustment 
is  required,  inverts  the  blocks,  and  uses  the  diagonal  elements  of  the 
inverses  to  calculate  the  error  degree  variances,  the  relative  error  per 
degree,  and  the  corresponding  accuracy  of  the  geoid  undulation;  the  second 
version  of  the  main  program  does  not  create  the  normals  but  reads  them 
from  tape  (unit  10,  once  more)  and  then  proceeds  as  before.  The  reason 
for  having  two  versions  .is  that,  with  the  first,  one  creates  the  normals 
and  then  carries  out  the  error  analysis  using,  say,  least  squares  adjustment 
theory;  as  the  normals  are  stored  without  C*1  ,  one  can  then  use  the 
second  program  to  add  C"1  to  the  stored  matrix  and  then  carry  out  the 
analysis  according  to  collocation  without  having  to  recompute  unnecessarily 
(AT  D"1  A)  ,  which  is  the  most  time-consuming  part  of  the  analysis. 

(a)  Full  Version 


This  program  both  creates  and  inverts  the  diagonal  blocks  of  the 
normal  required  by  the  analysis.  It  call  subroutines  ONEREV,  MATV ,  MODEL, 
NVAR,  and  WRIT  ; which  are  1 isted  in  this  Appendix;  it  also  calls  the  system- 
provided  subroutines  ERRSET,  SCL0K1,  and  RCL0K1,  and  the  subroutines  GGNOR 
(from  the  single  precision  library  produced  by  the  International  Mathematical 
and  Statistical  Libraries  Inc.  (IMSL)  company),  and  DSINV  from  the  double 
precision  library  of  the  IBM  System/360  Scientific  Subroutine  Package. 

The  program  contains  a  device  against  a  possible  underestimation 
of  the  total  execution  time  requested  in  the  JOB  card.  The  running  time 
is  checked  inside  the  loop  where  the  diagonal  blocks  are  created,  so  that 
if  this  time  is  a  few  seconds  below  the  assigned  running  time  before  a 
new  block  is  to  be  calculated,  the  run  terminates  in  an  orderly  way.  The 
blocks  created  so  far  are  left  safely  in  unit  10,  and  the  remaining  ones 
can  be  created  in, a  later  run  where  the  minimum  order  of  the  blocks  is 
set  to  MMIN  =  (order  of  last  (even,  odd)  pair  of  blocks  completed  previously 
+  1).  The  time  for  eventual  premature  exit  is  the  value  of  parameter 
YLIM  in  seconds,  and  should  be  smaller  than  the  time  declared  in  the  JOB 
card  by  at  least  10  seconds  +  compile  time.  Blocks  created  in  additional 
runs  should  be  stored  in  new  tapes,  which  can  be  combined  with  the  initial 
one  by  means  of  a  simple  FORTRAN  program  to  produce  a  tape  with  all 
blocks  in  the  required  sequence  (increasing  order).  The  running  time  is 
checked  with  the  help  of  SCL0K1  (which  sets  the  time  counter  to  zero  before 
the  main  do  loop)  and  RCL0K1  that  "reads"  the  time  counter  at  every  turn 
in  the  loop. 
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All  calculations  are  carried  out  in  double  precision,  with  the  excep¬ 
tion  of  those  related  to  subroutine  GGNOR,  that  take  place  in  single  pre¬ 
cision.  For  this  reason  all  variables  are  declared  REAL  8  ,  except  array 
TP  which  is  REAL  4  . 

The  only  COMMON  statement  is  "P",  through  which  the  value  of  the 
total  geoid  undulation  power  POWNMX  is  communicated  to  the  program  from 
subroutine  NVAR,  which  defines  the  values  of  the  degree  variances  according 
to  the  model  explained  in  section  3. 

The  various  arrays  are  somewhat  overdimensioned  for  the  needs  of 
the  actual  analysis,  when  N  =  331  .  The  maximum  degree  N  is  defined 
by  the  parameter  NMAX  which  must  be  an  odd  number,  for  reasons  given 
in  the  description  of  subroutine  ONEREV  (if  the  value  declared  happens 
to  be  even,  by  mistake,  the  program  adds  1  to  make  it  odd).  The  small 
arrays  (dimensioned  500)  should  be  at  least  of  dimension  NMAXP  =  NMAX 
+  1  .  As  for  the  large  arrays,  their  minimum  dimensions  should  be 

GPNMS  :  (NMAX  -  MMIN  +  2)  x  (NMAXP/2+1) 

AE  -g-  x  (NMAX  +  IV 

AO  : 

GPNMS  contains  the  Fourier  coefficients  of  the  columns  of  the  A  matrix, 
calculated  by  subroutine  ONEREV,  and  AE,  AO  contain  the  elements  of  the 
normal  matrix  corresponding  to  a  given  order  m  and  (n-m)  even  and 
odd,  respectively.  Storage  is  in  "upper-diagonal  form"  as  required  by 
the  inversion  subroutine  DSINV. 

Subroutine  ERRSET  is  used  to  suppress  unwanted  error  messages  due 
to  underflows  in  the  calculations  carried  out  by  ONEREV. 

The  run  parameters,  including  most  of  the  mission  parameters,  are 
assigned  values  between  statements  15  and  31.  Two  other  parameters,  the 
acceleration  of  gravity  at  ground  level  and  the  mean  Earth  radius,  are 
assigned  values  in  statements  50  and  51.  The  meanings  of  the  symbolic 
names  given  to  the  parameters  are  explained  in  the  comments  of  this  section 
of  the  program.  Parameter  IAC,  for  example,  tells  the  program  whether 
a  least  squares  adjustment  or  a  least  squares  collocation  error  analysis 
is  requested.  The  main  parameter  values  are  listed  by  the  program  (statements 
36  and  37  if  pTotter  output,  as  well  as  printer  output,  is  desired)  at 
the  beginning  of  the  run  and,  together  with  some  other  important  values, 
saved  at  the  start  of  the  file  created  in  unit  10,  where  the  diagonal 
blocks  are  also  stored.  The  length  of  day  is  upposed  to  be  24  x  3600 
seconds  exactly. 

Certain  loaders  initialize  all  array  areas  to  zero  before  execution, 
as  in  the  case  of  the  loader  used  to  run  this  program.  Others  set  arrays 
and  registers  not  declared  in  the  program  to  "indeterminate",  or  else 
the  values  left  in  core  by  a  previous  job  remain  unchanged.  In  such  cases 
the  program  does  not  work  unless  additional  DATA  statements  are  written 
to  give  all  array  elements  an  initial  value  of  zero. 
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The  run  begins  with  the  computation  of  sin  (4J-),  cos^),  and  other 
quantities  to  be  explained  later,  as  well  as  the  quantities  p  sinp^r 
cos  ««■  and  cosp-£  sin-j  stored  in  arrays  SPC  and  CPS,  respectively 
(statements  68  to 72).  this  numbers  are  used  later  to  form  the  Fourier 
coefficients  of  the  columns  of  A  according  to  (2.13).  Subroutines  NVAR 
and  MODEL  are  called  to  set  up  the  array  of  degree  variances  of  the  geopo¬ 
tential,  DVAR.  The  inverse  of  the  variances  of  the  scaled  coefficients 
(equations  (2.15,a-b))  are  stored  in  array  CM1.  The  variances  produced 
by  NVAR  correspond  to  n>  NMOD  and  are  in  the  form  of  gravity  anomaly 
variances.  They  are  converted  to  the  inverse  of  geopotential  variances 
in  statement  81.  All  information  relative  to  C_1  is  stored  in  CM1  and 
in  GMR. 

The  main  loop  begins  at  statement  91.  Statements  92  through  95  check 
that  the  running  time  allocated  in  the  JOB  card  is  not  exceeded.  The 
elapsed  time  from  the  beginning  of  the  main  loop  is  printed  at  the  beginning 
of  a  new  pass  through  the  loop  (statement  93).  Each  pass  creates  the 
two  diagonal  blocks,  one  for  n  even,  the  other  for  n  odd,  corresponding 
to  order  m  =  M  .  There  are  four  non-zero  blocks  for  each  m  ,  but  those 
corresponding  to  "sine"  terms  are  identical  to  those  for  "cosine"  terms 
(expression  (2.51)  is  independent  of  a  ),  so  only  one  pair  of  blocks 
is  needed.  Statements  97  through  105  calculate  the  factors 


FF<P>  t'1'  «»«*,— i*a)  <?£♦»£>- 

+  (1-  cos((pw-  mfi)Aa)(p-^-  -  ■J-)-*  ] 

U)0  W0 

(see  expression  (2.51))  where  WO  is  the  fundamental  angular  frequency 
for  the  whole  mission  (period  TO  =  Ndaysx24x3600) .  The  "apm"  calculated 
with  ONEREV  are  too  large  by  a  factor  of  2(N+1)  ,  and  this  fact  is  accounted 
for  in  the  denominator  above. 

After  ONEREV  has  been  called,  it  returns  the  "Fourier  coefficients" 


a™  *  hp1"  [  P  sin  P-|-  cos  +  (n+1)  cosp  sin 

where  the  hpm  correspond  to  all  Cnm^  ^  with  m  =  M  ,  n  z  M  ,  in 
array  GPNMS.  These  coefficients  are  then  multiplied  by  each  other  and 
scaled  by  the  factor  FF(P)  in  statements  107  through  137,  to  form  the 
element  g^nm  AIJ  according  to  (2.51).  If  a  collocation  analysis 
is  requested  and  AIJ  corresponds  to  an  element  on  the  diagonal,  the  cor¬ 
responding  term  in  C-1  is  added  to  AIJ  in  statement  136.  Between  138 
ana  154  the  calculated  elements  are  stored  in  AE  and  AO,  the  arrays  containing 
the  blocks  for  (n-m)  even,  and  (n-m)  odd,  respectively,  in  upper  diagonal 
form.  The  vectors  PIVE  and  PIVO  are  formed  with  the  square  roots  of  the 
diagonal  elements.  From  155  to  157  the  two  blocks  and  relevant  information 
regarding  their  size  are  stored  in  unit  10,  for  further  possible  use. 

From  158  to  165  the  two  blocks  are  "pivoted"  according  to 


rP  reven,  _ 
Va,  ‘odd  »  ~ 


p-i 


r  reveni 
m,a,  ‘odd  f 
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P'1 


(B.l) 


where  P  is  a  matrix  of  "pivots",  i.e.,  the  square  roots  of  the  diagonal 
elements  of  Gm.a.iodd  and  P I  VO  are  used  for  pivoting  AE  and 

AO  ,  respective fy. .  The  "pivoted"  blocks  are  inverted  in  statements  173 
and  174  (there  is  no  "AO"  block  for  M  =  N ) .  Notice  that  the  program  has 
been  written  in  such  a  way  that  AE  and  AO  correspond  to  even  and  to  odd 
n-m  ,  not  n  .  As  m  is  fixed  during  each  pass  of  the  main  loop,  the 
result  depends  on  n  ,  so  AE  sometimes  contains  "even",  and  sometines 
"odd"  elements,  and  the  reciprocal  is  true  of  AO,  depending  on  the  value 
of  m  . 

The  accuracy  of  the  inverse  is  tested  by  the  inversion  subroutine 
DSINV  according  to  the  tolerance  limit  defined  in  the  main  program  (statement 
172).  If  DSINV  returns  a  value  of  zero  in  register  IER1,  the  inversion 
can  be  regarded  as  successful  (free  of  serious  numerical  instability). 

In  addition  to  this  test  by  DSINV,  a  further  check  has  been  written  in 
the  segment  from  statement  167  to  194.  The  idea  is  to  calculate 


«t  *  t  -kg£  r1  Sn)]t-t-f 

-  -  "i,«,  'odd  (computed)  °dd  '  '  - 


and  then 


[6tT  fit  (tT  t)_1  ]* 


The  exponent,  in  floating  point  notation,  of  e  is,  approximately,  the 
number  of  significant  figures  within  which  t  and  t1  agree.  Unfortunately,, 
because  188  and  192  are  not  coded  correctly ,*the  values  of  e  printed 
at  statement  193  are  not  useful.  The  elements  of  t  are  random  numbers 
created  by  the  IMSL  subroutine  GGNOR. 

From  196  to  the  end  of  the  main  loop  the  inverted  normal  blocks  are 
obtained  by  "de-pivoting" : 

g;\  {!!!")  =  p'1  (g!  _  f!^ni  r1  p-1  (b.2) 


and  the  variance  of  the  error  per  degree  is  totalized  in  array  RMS  as 
follows 


RMS(„)  =  RHS(n)  *  m  (f )!"  x  ( 


1  if  m  =  0 

2  otherwise 


(o2Cnmnm  =  o2£nmnm)  so>  at  the  end  of  the  ma1n  loop.RMS(n)  =  -jT*(|-)'n'72en  , 
where  o2en  is  the  n  degree  variance  of  the  errors,  and  the  rest  is 
the  scaling  factor  squared  (expressions  ( 2 . 15 ,a -b ) ) .  From  statement  219 
to  the  end  the  error  in  the  undulation  up  to  degree  n  and  this  error 
plus  the  truncation  error  above  n  are  both  calculated  and  stored  in 
arrays  PPB  and  PPTR,  respectively.  Array  PERC  receives  the  formal  percen¬ 
tage  error  per  degree  (expression  (3.2)  multiplied  by  100).  While  "debug¬ 
ging"  the  program  it  was  found  useful  to  monitor  the  values  of  some  of 
the  coefficients'  accuracies  as  they  were  being  obtained.  This  feature 
was  left  in  the  program,  where,  for  the  reason  given  in  paragraph  3.4, 
all  the  o2e$^nm  for  n  s  40  are  still  printed  out  (statement  210). 

If  the  inversion  of  AE  or  of  AO  fails,  subroutine  DSINV  returns  an  explanatory 


code  in  IER1  and  IER2,  which  should  be  different  from  0  and  is  also  printed 
out.  A  failure  to  invert  is  most  likely  due  to  a  set-up  error;  perhaps 
to  an  improperly  dimensioned  array;  to  a  loader  that  does  not  preset  all 
undefined  values  to  0  before  execution;  or  to  improperly  punched  cards. 

An  example  of  all  the  messages  printed  during  a  normal  run  can  be  seen 
in  paragraph  B.4. 

The  program  prints  out  all  results  in  unit  III  (a  plotter,  for  instance) 
degree  by  degree  (statement  238)  in  blocks  of  fifty  lines,  and  at  the 
end  prints  a  summary  in  increments  of  10  degrees.  The  results  are  also 
punched  out  (statements  219  and  220). 

(b)  Reduced  Version 

This  version  of  the  main  program  does  not  compute  the  normals,  but 
reads  them  from  unit  10  (disk  or  tape),  where  they  have  been  stored  during 
apreviousrun  of  the  full  version.  It  also  reads  the  original  mission 
parameters,  with  which  the  normals  were  created,  as  the  first  record  in 
file  10  (stat.  22).  Some  of  those  parameters  can  be  changed  in  value 
by  re-scaling  the  normals  (paragraph  2.13).  The  new  parameters'  values 
can  be  declared  by  inserting  statements  between  lines  22  and  25  in  the  listing. 

This  version  does  not  call  subroutine  ONEREV,  and  reads  unit  10  from 
subroutine  RED.  Subroutine  WRIT  is  not  used.  All  the  other  subroutines 
called  in  the  full  version  are  also  used  here.  The  fact  that  there  is 
no  AO  matrix  block  when  M  =  NMAX  is  taken  into  account  (statements 
94  and  98).  Re-scaling  according  to  (2.78)  happens  between  99 
and  104.  If  no  change  in  parameters  is  desired,  FNSR  is  1  .  Expression 
(2.84)  for  least  squares  collocation  requires  knowledge  of  the  singular 
values  and  eigenvectors  of  the  normal  matrix,  instead  of  the  inverse. 

Lack  of  time  and  of  familiarity  with  the  decomposition  subroutines  available 
resulted  in  inversion  being  chosen  for  both  least  squares  adjustment  and 
for  collocation.  The  fastest  way  to  study  the  effect  of  parameter  value 
changes  on  the  results  with  collocation  is,  therefore,  to  run  this  reduced 
version  with  the  new  parameters.  With  NMAX  =  331,  this  requires  some 
15  minutes.  As  no  PIVE  or  PIVO  arrays  have  been  computed  so  far,  this 
is  done  now  in  statements  122-125,  and  then  pivoting  is  applied.  From 
there  on  things  are  done  much  as  in  the  full  version,  except  that  no  test 
for  numerical  stability  of  the  inversion  is  done  in  addition  to  that  of 
DSINV.  Results  are  printed  and  punched  as  in  the  full  version. 


B.2  Subroutine  ONEREV 

This  subroutine  computes  the  coefficients  apm  .teeded  to  fin 
elements  of  the  normal  matrix  according  to  expression  (2.51).  The  a. 
are  defined  by  expression  (2.13) 

apm  =  h£m  [(n  +  1)  cos  p  sin  +  p  sin  p  ^  cos  y  ] 

as  proportional  to  the  Fourier  coefficients  Tipm  of  the  extended  Legendre 
functions  Cnm^')  •  The  values  of  (n+1)  cos  p  ^  sin-*r  and  of  p  sinp| 
cos are  passed  on  to  the  subroutine  from  the  main  program  in  arrays  SPC 
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and  CPS*  respectively.  The  Cnm(4  )  are  sampled  at  equal  intervals 
A<f>*  =  i  (N  is  the  highest  degree  in  the  band,  so  the  highest  fre¬ 

quency  Temrin  any  Cnm  is  either  hfl'1’' cosN  $ '  or  h^m  sinN^‘)»  and 
then  analysed  with  subroutine  FFCSIN,  which  implements  a  mixed-radix  Fast 
Fourier  Transform  (FFT)  algorithm.  FFCSIN  belongs  to  the  International 
Mathematical  and  Statistical  Libraries  (IMSL)  Inc.'s  double  precision 
library.  It  returns  2{N+l)T)pm  instead  of  hpm  ,  but  this  is  taken  into 
account  in  the  main  program.  The  Lnm^  '  are  calculated  taking  advantaqe 
of  the  relationships  {2.2,  a-d),  so  only  the  values  in  the  interval 
0  <  <j>  <  are  absolutely  necessary.  These  values  are  the  same  as  those 
of  the  corresponding  Pnm  ,  which  are  obtained  with  subroutine  LEGFDN.  ^ 
There  are  (N+l)/2  points!^  where  the  Pnm  have  to  be  calculated  in  0  £$<-*-, 
so  the  space  required  by  GPNMS  is  (N  -  MMIN  +2)  x  ( (N  +  1/2). 

Here  MMIN  is  the  lowest  order  to  be  studied  (a  feature  of  the  main  program 
is  that  it  allows  the  study  of  coefficients  in  the  band  MMIN  <  m  s.  N  ; 
in  the  case  of  the  results  of  section  3,  MMIN  =  0).  The  successive  values 
of  the  Pnm  are  put  first  in  array  GPNMS.  Then  all  values  corresponding 
to  the  same  Pnm  are  moved  to  REV,  where  Cnm  is  determined  from  (2. 2, a-d), 
so  the  dimension  of  REV  js  2(N+1).  Then  FFCSIN  replaces  the  values  of 
Cnm  with  those  of  the  hp  ,  also  in  REV,  and  the  latter  are  moved  on 
to  GPNMS,  where  they  replace  the  original  Pnm  •  IWK  is  an  auxiliary 
array  needed  by  FFCSIN  (see  IMSL  Handbook).  All  aj3m  (multiplied  by 
2(N+1))  are  returned  to  the  main  program.  If  the  loader  used  does  not 
preset  all  undefined  registers  to  zero,  a  DO  loop  should  be  put  at  the 
beginning,  between  statements  5  and  6,  setting  to  zero  all  arrays  with 
dimensions  different  from  1  .  Alternatively,  depending  on  the 

compiler,  a  DATA  statement  to  the  same  effect  could  be  used. 

If  an  even  function  is  added  to  an  odd  function,  the  Fourier  coefficients 
of  the  sum  are  those  of  the  even  function  in  the  cosine  terms,  and  those 
of  the  odd  function  in  the  sine  terms.  This  simple  fact  is  exploited 
to  reduce  calculations  by  half,  taking  advantage  of  the  even  or  odd  nature 
of  the  Lnm  with  respect  to  O'  . 

The  a[)m  with  p=0  are  handled  separately  from  the  rest,  and  are 
returned  to  the  main  program  in  array  GMN  (statement  78).  Besides  FFCSIN 
and  LEGFDN,  no  other  subroutines  are  called. 


B.3  Subroutines  LEGFDN,  MODEL,  and  NVAR 


Subroutine  LEGFDN  can  calculate  both  the  values  of  all  normalized 
Legendre  function^  at  colatitude  THETA,  for  the  same  order  m  =  M  ,  up 
to  degree  n  =  NMAX  ,  and  also  their  derivatives.  The  functions  are  returned 
in  array  RLEG,  the  derivatives  in  DLEG ;  after  exit,RLMN  contains  a  redundant 
set  of  all  sectorial s  up  to  degree  n  =  m  .  All  arrays,  except  DRTS 
and  DIRT  should  have  the  dimension  NMAX1  (statement  5).  DRTS  and  DIRT 
should  have  twice  that  size.  IR  is  a  register  that  should  be  set  to  zero 
before  the  first  call  to  this  subroutine,  in  the  main  program.  IFLAG 
tells  the  subroutine  whether  only  the  Legendre  functions  or  these  and 
their  derivatives  are  required  (only  the  Pnm  were  needed  for  the  error 
analysis,  so  IFLAG  was  set  to  l).  The  Pn«  and  their  derivatives  are 
calculated  using  recursive  formulas  given  in  Colombo  ((1981),  paragraphs 
(1.10)  and  (4.4)). 

^(N  +  l)/2  must  be  integer,  so  N  »  NMAX  must  be  an  odd  number. 
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Subroutine  MODEL  sets  the  values  of  the  first  NMOD  components  of  arrays 
DVAR  adn  DVER  to  the  values  of  the  degree  variances  of  the  errors  in  the 
potential  coefficients  of  the  reference  model  with  respect  to  which  the 
residual  line  of  sight  velocities  are  determined.  The  values  in  the  listing 
corespond  to  the  first  30  degrees  in  a  model  obtained  by  R.  H.  Rapp  at 
OSU  from  a  global  data  set  of  l°xl°mean  anomalies  by  numerical  quadratures. 
This  model  is,  in  fact,  complete  up  to  degree  and  order  180,  though  only 
the  accuracies  of  the  first  30  degrees  are  used  here. 

Subroutine  NVAR  initializes  array  DVAR  so  that  its  elements  are  the 
same  as  the  degree  variances  for  the  gravity  anomaly  implied  by  R.  Rapp's 
coefficients  mentioned  above,  up  to  n  =  100  .  Above  n  =  100  ,  the  variances 
are  those  obtained  from  a  two-term  model  for  the  gravity  anomaly  spectrum, 
also  the  work  of  R.  Rapp  (1979b).  The  main  program  calls  NVAR  first  and 
MOOEL  afterwards,  so  the  first  NMOD  degree  variances  are,  finally,  those 
in  MODEL  (error  variances),  which  correspond  to  dimensionless  potential 
coefficients.  The  remainder  comes  from  NVAR,  so  they  must  be  converted 
to  dimensionless  potential  from  gravity  anomaly  variances,  a  step  that 
occurs  in  statements  80  and  81  of  the  main  program.  In  NVAR,  statements 
120  and  122  add  all  potential  degree  variances  between  n  *  NMAX  and  n 
=  2000  ,  and  return  this  sum  in  POWNMX  through  COMMON/P/  to  the  main 
program. 

The  same  caution  regarding  the  initialization  of  undeclared  arrays 
and  variables  to  zero  that  was  made  for  the  main  program  and  for  0NEREV 
apply  to  LEGFDN  and  to  MODEL  (array  DVAR)  as  well. 


B.4  Sample  Output 

After  the  listings  of  the  various  routines  described  previously, 
the  reader  will  find  a  sample  of  the  printed  output  created  by  either 
version  of  the  main  program  (they  are  the  same).  This  listing  should 
help  whoever  wants  to  use  these  programs  to  check  that  his  own  punched 
version  works  properly.  The  program  should  also  punch  out  some  cards 
with  the  results  (statements  238  and  239).  Because  of  a  minor  error  in 
the  program,  now  corrected,  the  standard  deviation  of  the  data  is  listed 
as  10*6  m  s*1,  instead  of  as  /Z  x  10-6  m  s-1,  which  is  the  actual  value 
corresponding  to  the  results  printed  in  the  sample.  The  parameters  used, 
with  the  exception  of  o  ,  are  as  listed.  To  enter  them  into  the  main 
program,  see  comments  at  the  beginning  of  either  version. 

The  first  page  of  the  output  contains  a  listing  of  the  mission  parameter 
values  chosen.  The  "TIME  BEFORE  M  ”  statements  give  the  time  in  seconds 
at  the  beginning  of  a  new  pass  through  the  main  loop,  and  they  are  printed 
just  before  the  time  check  in  statement  95.  the  "ACCURACY  OF  INVERSION" 
lines  should  indicate  the  stability  of  the  inversion  of  the  blocks  of 
order  M  ,  but  they  are  useless  because  of  the  coding  error  mentioned 
in  paragraph  B.l.  The  other  numbers  are  scaled  variances 
o2enm«m  of  the  coefficients  up  to  degree  and  order  40  (a=0)  .  If  the 
numerical  inversion  of  either  block  of  order  M  fails  the  stability  test 
in  subroutine  DSINV,  a  line  is  printed  saying  "AT  ORDER  m  IER1  =  x 
1ER2  *  y",  where  x  and  y  are  two  integers  whose  value  should  be  interpreted 
according  to  instructions  in  the  Handbook  of  the  SSP  library.  At  the  end 
of  the  run,  the  various  accuracies  are  listed  in  unit  IU,  first  degree 
by  degree  and  then,  in  a  final  summary  page,  every  10  degrees. 
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Appendix  C  :  Detailed  Listings  Degree  by  Degree 
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